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

Question

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. — Henrik Schumacher, Dec 24 '19 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. — m_goldberg, Dec 24 '19 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? — m_goldberg, Dec 24 '19 at 16:03
You could apply FindRoot with randomly chosen starting points. You might describe that as Monte Carlo. — mikado, Dec 24 '19 at 19:10
@m_goldberg Thanks. I can only find the reference in Chinese, https://www.ixueshu.com/document/dffe3e67c253298a318947a18e7f9386.html. If you are interested in it and can not open the URL, I can email to you. — keanhy14, Dec 25 '19 at 01:51
I'm afraid that a reference in Chinese is beyond me. I am totally ignorant of Chinese. — m_goldberg, Dec 25 '19 at 02:15
@All, this paper presents a method, and the first example is a univariate nonlinear equation: https://www.sciencedirect.com/science/article/pii/0771050X80900224 — Michael E2, Dec 25 '19 at 06:13

Xminer · Answer 1 · 2019-12-28T08:46:52.300

2

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}]

Mathematica graphics

MAP estimating

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

{x -> 1.10409}

edited Dec 28 '19 at 08:46

answered Dec 26 '19 at 13:14

Xminer

1
7
15

Thanks! Is it convenient to provide some materials (book or article) to understand your code? – keanhy14 Dec 27 '19 at 01:46
I'd like to recommend course slides(generative models,1~4) of "Machine Learning for Engineers" – Xminer Dec 27 '19 at 03: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]? – Michael E2 Dec 27 '19 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. – Xminer Dec 28 '19 at 08:49

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

1 Answers1