restart

with(plots):
with(Statistics):

# 1st example
#
# It is a simple 1D example for the uniform measure seems the natural one

f := x -> exp(x):

N := 15000:
a := 0:
b := 1:

X := Statistics:-Sample(Uniform(a, b), N):
(b-a)*Mean(f~(X));

evalf(Int(f(x), x=0..1));

# 2nd example (yours)
#
# It is a 2D example in (what seems) polar coordinates.
# A more complicated example that the first one wich forst implies defining what is the good probability measure to use.


f := (rho, phi) -> exp(rho*cos(phi));


T3D    := plottools:-transform((x, y, z) -> [x*cos(y), x*sin(y), z]):
cart3D := plot3d(exp(rho*cos(phi)), rho = 0 .. 1, phi = 0 .. 2*Pi, style=surface):

DocumentTools:-Tabulate([cart3D , T3D(cart3D)], width=60)

# Do we have to consider this "2D cartesian" cloud of points (left figure) which transforms
# in this "2D polar" cloud (right figure) ?

T := plottools:-transform((x, y) -> [x*cos(y), x*sin(y)]):

N := 1000:
X := Sample(Uniform(0, 1), N):
Y := Sample(Uniform(0, 2*Pi), N):

cart2D  := ScatterPlot(X, Y, labels=[typeset(rho), typeset(phi)]):

DocumentTools:-Tabulate([cart2D , display(T(cart2D), polarplot(1))], width=60)

# Or do we have to construct a "2D polar" cloud with an ad hoc uniformity spreading?
# for instance a distribution of points that preserve the mass (Pi) of the unit disc?


PHI := Sample(Uniform(0, 2*Pi), N):
RHO := sqrt~(Sample(Uniform(0, 1), N)):

X := RHO *~ cos~(PHI):
Y := RHO *~ sin~(PHI):
polar2D := ScatterPlot(X, Y):
display(polar2D, polarplot(1))
 

# Now let's integrate the finction f in polar coordinates
#
# Found by Maple's implicit integration method

evalf( Int(exp(rho*cos(phi))*rho, rho = 0 .. 1, phi = 0 .. 2*Pi) );

# Found by Maple with method=_MonteCarlo or Cuba variants, see here evalf/Int
# (cancelled after 30 seconds, but it's likely the result would have been around 5.4 because
# the methode operates on hyperrectangles)

# evalf( Int(exp(rho*cos(phi))*rho, rho = 0 .. 1, phi = 0 .. 2*Pi, method = _CubaDivonne) );

Warning,  computation interrupted

# 1st 2D measure based form a Bi-Uniform cartesian measure... which is probaly the wrong one

N   := 15000:
X   := Sample(Uniform(0, 1), N):
Y   := Sample(Uniform(0, 2*Pi), N):
RHO := sqrt~(X^~2 + Y^~2):
PHI := arctan~(Y/~X):

area := Pi:

FirstValue := area*Mean(f~(RHO, PHI));

# 2nd 2D measure based form a a uniform sampling of the unit disc... which is probaly the good one
# (at least the more natural) unless otherwise stated

N   := 15000:
PHI := Sample(Uniform(0, 2*Pi), N):
RHO := sqrt~(Sample(Uniform(0, 1), N)):

area       := Pi:
FirstValue := area*Mean(f~(RHO, PHI));

# As you probably know the speed of convergence of a Monte Carlo (MC) method is N^(-1/2).
# More precisely the variance of the error between the MC summation and the exact integration
# behaves as N^(-1).

MC := proc(R, N)
  local E, r, PHI, RHO:
  E := Vector(R):
  for r from 1 to R do
    PHI  := Sample(Uniform(0, 2*Pi), N):
    RHO  := sqrt~(Sample(Uniform(0, 1), N)):
    E[r] := area*Mean(f~(RHO, PHI));
  end do:

  Variance(E);
end proc:

CodeTools:-Usage( MC(1000, 100) );

memory used=46.60MiB, alloc change=0 bytes, cpu time=967.00ms, real time=3.42s, gc time=104.14ms

CodeTools:-Usage( MC(1000, 1000) );

memory used=184.28MiB, alloc change=256.00MiB, cpu time=4.41s, real time=13.44s, gc time=1.99s

CodeTools:-Usage( MC(1000, 10000) );

memory used=1.52GiB, alloc change=32.00MiB, cpu time=29.99s, real time=27.07s, gc time=13.94s

# There exist MC like methods with faster (asymptotic) convergence rate.
# For instance Quasi Monte Carlo (QMC) methods, but they "work" only for integration on hyperrectangles
# (not your case).
# Latin Hypercude Sampling (LHS) methods have a better behaviour than MC methods for they
# produce a more homogenous distribution of samples than the latter.
#
# It has been prooved that the integration error is of the form V(f)*D(X) (Koksma–Hlawka inequality)
# https://en.wikipedia.org/wiki/Low-discrepancy_sequence
# where V(f) is a term (Hardy-Krause's total variation) which depens only on the function f to integrate,
# and D(X) a term (called Discrepancy) that depends only on the summation sample X.
#
# Reducing the error implies reducing this discrepancy.
# QMC, Minimax, Maximin, CVT (Centroïdal Voronoï Tesselation)... and many others, are methods with a
# smaller asymptotic discrepancy (for a same number of points) than crude MC method.
#
# The term asymptotic is extremely important because all these methods, but MC, contain a
# multiplicative coefficient which grows exponentially with the multiplicity of the integral,
# making their theoritical advantage more questionable in high dimensions.

restart;

with(Statistics):

X := [-0.012, -0.010, -0.004, -0.002, -0.001, -0.0001, 0.0001, 0.001, 0.002, 0.004, 0.010, 0.012]:
Y := [-0.695, -0.7, -0.74, -0.825, -0.95, -1.0, 1.0, 0.95, 0.825, 0.74, 0.7, 0.695]:

MM := x -> (a+b*ln(c*abs(x)))*signum(x)

Target := add((Y -~ MM~(X))^~2):
opt_M  := Optimization:-Minimize(Target, {c >= 0});

MF := unapply(eval(MM(x), opt_M[2]), x)

plots:-display(
        plot(MF(x), x=min(X)..max(X), thickness=3),
        plots:-pointplot(X,Y,symbol=solidcircle,symbolsize=12,color=blue),
        view=-2..2, size=[500,350]
);

restart;

with(Statistics):

X := [-0.012, -0.010, -0.004, -0.002, -0.001, -0.0001, 0.0001, 0.001, 0.002, 0.004, 0.010, 0.012]:
Y := [-0.695, -0.7, -0.74, -0.825, -0.95, -1.0, 1.0, 0.95, 0.825, 0.74, 0.7, 0.695]:

# from acer's reply

Xn,Xp := selectremove(`<=`,X,0):
T := table([seq(X[i]=Y[i],i=1..nops(X))]):
Yp,Yn := [seq(T[x], x=Xp)], [seq(T[x], x=Xn)]:

Mn := convert([Xn, Yn], Matrix)^+:
Mp := convert([Xp, Yp], Matrix)^+:

# Observe the symmetry of the data


Mn, Mp[sort(Mp[..,2], output= permutation)];

# Viviane's model : y=c+exp(-b*x)


Fp := NonlinearFit(c+exp(-b*x), Mp, x):
Fp := unapply(Fp, x);

# Fn respects the symmetry of the data

Fn := x -> -Fp(-x)

# Residual Sum of Squares

RSS_p := add((Yp -~ Fp~(Xp))^~2);
RSS_n := add((Yn -~ Fn~(Xn))^~2);

P := plots:-pointplot(X,Y,symbol=solidcircle,symbolsize=12,color=blue):
plots:-display(
        plot(Fn(x), x=min(Xn)..max(Xn), thickness=3),
        plot(Fp(x), x=min(Xp)..max(Xp), thickness=3),
        P,
        view=-2..2, size=[500,350]
);

# Acer's model y = a+b*ln(abs(c*x))

AM := x -> a+b*ln(c*abs(x));

AM_p := x -> a+b*ln(c*x);
AM_n := x -> a+b*ln(-c*x);
 

Mn[sort(Mn[..,2], output= permutation)], Mp

tol := optimalitytolerance=1e-10:

Target_p := add((Yp -~ AM_p~(Xp))^~2):
opt_p    := Optimization:-Minimize(Target_p, tol);

Target_n := add((Yn -~ AM_n~(Xn))^~2):
opt_n    := Optimization:-Minimize(Target_n, tol);




Target_p0 := add((Yp -~ AM~(Xp))^~2):
opt_p0    := Optimization:-Minimize(Target_p0, tol);

Target_n0 := add((Yn -~ AM~(Xn))^~2):
opt_n0    := Optimization:-Minimize(Target_n0, tol);

# Lack of antisymmetry between opt_n and opt_p (or opt_n0 and opt_p0) come probably from different
# initial points.
# Here are the solutions obtained for initial points that are "antisymmetric")

opt_p  := Optimization:-Minimize(Target_p , tol, initialpoint = {a= 0.4, b=-0.1, c=1});
opt_n  := Optimization:-Minimize(Target_n , tol, initialpoint = {a=-0.4, b= 0.1, c=1});
opt_p0 := Optimization:-Minimize(Target_p0, tol, initialpoint = {a= 0.4, b=-0.1, c=1});
opt_n0 := Optimization:-Minimize(Target_n0, tol, initialpoint = {a=-0.4, b= 0.1, c=1});

restart

a := 5:
b := 7:
k := 9:
A := [a, 0]:
B := [0, b]:

#A et B fixes

P := [t, 0]:
Q := [0, k/t]:

#P et Q 2 points mobiles

cir := -a*x-b*y+x^2+y^2 = 0;
sol := solve(subs(y = 5, cir), x);
cen := [solve(diff(cir, x)), solve(diff(cir, y))];
x0  := sol[1]:
y0  := 5:
M   := [x0, y0]:
R   := sqrt(cen[1]^2+cen[2]^2);


EQ(M, cen); #... EQ is undefined...

#... thus solve returns nothing

solve(EQ(M, cen), y)

beta := arctan(diff(solve(EQ(M, cen), y), x)); #... and diff(nothing, x) generates an error

Error, invalid input: diff expects 2 or more arguments, but received 1

# Recherche des valeurs de t pour que les 2 droites soient perpendiculaires

eq := t^2*(y0-b)+t*(a*b-a*y0+b*x0-k)-x0*(a*b-k) = 0;
sol := solve(eq, t);
t := sol[1];
tp := sol[2];
P1 := [t, 0];
Q1 := [0, k/t];
PQ1 := simplify(x*(-a*b+b*t+k)+y*t*(t-a)-t*(-a*b+b*t+k)) = 0:

#1ere tangente

PQ2 := simplify(x*(-a*b+b*tp+k)+y*tp*(tp-a)-tp*(-a*b+b*tp+k)) = 0:

#2ième tangente P2 := [tp, 0]; Q2 := [0, k/tp];

CIR := implicitplot(cir, x = -4 .. 8, y = -4 .. 12, color = red);
 

Error, invalid input: diff expects 2 or more arguments, but received 1

Warning, solving for expressions other than names or functions is not recommended.

Error, (in solve) a constant is invalid as a variable, 5/2+(1/2)*65^(1/2)

Error, invalid subscript selector

Error, invalid subscript selector

Fig := proc (alpha)
  local Dr1, DR1, Dr2, DR2, N, u0, v0, Po, t, tp, sol;
  global a, b, k, cen, R;
  u0  := cen[1]+R*cos(alpha);
  v0  := cen[2]+R*sin(alpha);
  N   := [u0, v0];
  sol := solve(t^2*(v0-b)+t*(b*u0-a*v0+a*b-k)-u0*(a*b-k) = 0, t);
  t   := sol[1];
  tp  := sol[2];
  Dr1 := simplify(x*(-a*b+b*t+k)+y*t*(t-a)-t*(-a*b+b*t+k)) = 0;
  DR1 := implicitplot(Dr1, x = -4 .. 8, y = -4 .. 12, color = brown);
  Dr2 := simplify(x*(-a*b+b*tp+k)+y*tp*(tp-a)-tp*(-a*b+b*tp+k)) = 0;
  DR2 := implicitplot(Dr2, x = -4 .. 8, y = -4 .. 12, color = pink);
  Po  := pointplot([N[]], symbol = solidcircle, color = [black], symbolsize = 8);
  plots:-display([Po, DR1, DR2])
end proc;
 

DrPQ1 := implicitplot(PQ1, x = -4 .. 22, y = -4 .. 12, color = blue);
DrPQ2 := implicitplot(PQ2, x = -4 .. 22, y = -4 .. 12, color = blue);
Points := pointplot([A[], B[], M[], P1[], P2[], Q1[], Q2[], cen[]], symbol = solidcircle, color = [green], symbolsize = 10);
T := plots:-textplot(
          [[A[], "A"], [B[], "B"], [M[], "M"], [P1[], "P1"], [P2[], "P2"], [Q1[], "Q1"], [Q2[], "Q2"], [cen[], "cen"]],
          font = [times, 10], align = {below, left}
     );
n := 19;

plots:-display([seq(Fig(2*i*Pi/n), i = 0 .. n), Fig(beta), CIR, DrPQ1, DrPQ2, Points, T], scaling = constrained, size = [500, 500])

Error, invalid input: diff expects 2 or more arguments, but received 1

Error, (in plots:-textplot) improper op or subscript selector

Error, (in plots:-display) cannot make plot structure from object with name pointplot