9
$\begingroup$

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).

Improve this question
$\endgroup$
8
  • $\begingroup$ 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. $\endgroup$
    – Henrik Schumacher
    Dec 24, 2019 at 15:05
  • 2
    $\begingroup$ 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. $\endgroup$
    – m_goldberg
    Dec 24, 2019 at 16:00
  • $\begingroup$ Can you supply a reference to a text book or article on using Monte Carlo methods to solve a single non-linear equation? $\endgroup$
    – m_goldberg
    Dec 24, 2019 at 16:03
  • 1
    $\begingroup$ You could apply FindRoot with randomly chosen starting points. You might describe that as Monte Carlo. $\endgroup$
    – mikado
    Dec 24, 2019 at 19:10
  • 3
    $\begingroup$ @All, this paper presents a method, and the first example is a univariate nonlinear equation: sciencedirect.com/science/article/pii/0771050X80900224 $\endgroup$
    – Michael E2
    Dec 25, 2019 at 6:13

1 Answer 1

Reset to default
2
$\begingroup$

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}
Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks! Is it convenient to provide some materials (book or article) to understand your code? $\endgroup$
    – keanhy14
    Dec 27, 2019 at 1:46
  • $\begingroup$ I'd like to recommend course slides(generative models,1~4) of "Machine Learning for Engineers" $\endgroup$
    – Xminer
    Dec 27, 2019 at 3:01
  • $\begingroup$ 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]? $\endgroup$
    – Michael E2
    Dec 27, 2019 at 15:04
  • $\begingroup$ @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. $\endgroup$
    – Xminer
    Dec 28, 2019 at 8:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.