# How can I solve a nonlinear equation using a Monte Carlo method?

I know FindRoot can solve nonlinear equations, but now I'm interested in solving such an equation using a Monte Carlo method. Can somebody show how that could be done for the equation below?

Exp[-x^3] - Tan[x] + 800 == 0


For convinience：The root is restrained in the interval (0,π/2).

• I do not know of any built-in routine that can do this. I guess you have to look into the literature, dig out one of the algorithms there, and implement it yourself. Dec 24, 2019 at 15:05
• I have heard of Monte Carlo methods being applied to systems of linear equations but never as a way to solve a single non-linear equation. I think that it is highly unlikely that there is anything built into to Mathematica that will apply such methods to a non-linear equation. Dec 24, 2019 at 16:00
• Can you supply a reference to a text book or article on using Monte Carlo methods to solve a single non-linear equation? Dec 24, 2019 at 16:03
• You could apply FindRoot with randomly chosen starting points. You might describe that as Monte Carlo. Dec 24, 2019 at 19:10
• @All, this paper presents a method, and the first example is a univariate nonlinear equation: sciencedirect.com/science/article/pii/0771050X80900224 Dec 25, 2019 at 6:13

Isn't basically just setting the likelihood and sampling x-values?
Let's define f[x]

f[x_] := Exp[-x^3] - Tan[x] + 800;


(Just Solve for comparing with monte-carlo)

NMinimize[{Abs[f[x]], 0 <= x <= Pi/2}, x]


{7.34717*10^-6, {x -> 1.56955}}

Think Distribution of $$y=f[x]$$.
Delta Method can be used.

joint[x_, y_, sigma_] :=
PDF[NormalDistribution[f[x], Evaluate@D[f[x], x]*sigma], y];


because here we have $$y=0(f[x]=0)$$,we can obtain the likelihood $$L[x]$$

L[x_?NumericQ, sigma_?NumericQ] := Evaluate@joint[x, 0, sigma];


Let's sampling x-values with mathematica-mcmc

Import["https://raw.githubusercontent.com/joshburkart/mathematica-\
mcmc/master/mcmc.m"]

SeedRandom[1234];
mcmc = MCMC[
L[x, sigma], {{x, 0.01, 0.01, Range[0, Pi/2, 0.01]}, {sigma, 100,
100, Range[100, 100000, 100]}}, 10000000]


obtain distribution of x-values

dst = mcmc["ParameterRun"][[;; , 1]] // SmoothKernelDistribution;
Plot[{Evaluate@PDF[dst, x]}, {x, 0, 2}]


MAP estimating

Last@NMaximize[Evaluate@PDF[dst, x], x]

{x -> 1.10409}

• Thanks! Is it convenient to provide some materials (book or article) to understand your code? Dec 27, 2019 at 1:46
• I'd like to recommend course slides(generative models,1~4) of "Machine Learning for Engineers" Dec 27, 2019 at 3:01
• Isn't the problem with the accuracy explained by the large derivative at or near the minimum/root, which is next to the asymptote of f[x]? Dec 27, 2019 at 15:04
• @MichaelE2 At First,I felt that there was a problem with accuracy because the convergence destination differs depending on the given initial interval. However, it doesn't seem to make much difference for MAP estimation if given sample is enoguh. So I corrected the article. Dec 28, 2019 at 8:49