2

I have the following function 

f(x)∝|x| exp(-1/2 |x| )+1/(1+(x-40)^4 ),xϵR

I want to find out E(X) and E(X^3) through Simpson's method (numerical integration), Standard Monte Carlo approach, Acceptance-rejection sampling, Importance sampling, Metropolis-Hasting Algorithm, Gibbs sampling and then Bayesian model using MCMC (I have not decided yet).

How can I validate my results obtained from different methods? I have tried to solve E(X) mathematically but fail to find any close form. This function can be divided over different parts as

absolute(x)*double exponential density + another function utilizing higher power (4) of X in inverse form. Due to absolute (x) and range [-Inf, Inf] We always have to divide it over [-Inf, 0] and [0, Inf]. Through Integration by parts I was able to see first part as (absolute (x) + (x^2/2) over infinite range) + Integral of this part can't be found mathematically.

So I make use of the following code to get numerical integration result as 

Library(stats)
integrand <- function(x) {x*(abs(x)* exp(-0.5*abs(x))+(1/(1+(x-40)^4)))}
integrate(integrand, lower = -Inf, upper = Inf)

thus the result is E(X)= 88.85766 with absolute error < 0.004

The results which I obtain from these methods are not similar for instance

(i) Through Simpsons method I got E(X) = 0.3222642 and E(X^3)=677.0711.. 

simpson_v2 <- function(fun, a, b, n=100) {
    # numerical integral using Simpson's rule
    # assume a < b and n is an even positive integer
    if (a == -Inf & b == Inf) {
        f <- function(t) (fun((1-t)/t) + fun((t-1)/t))/t^2
        s <- simpson_v2(f, 0, 1, n)
    } else if (a == -Inf & b != Inf) {
        f <- function(t) fun(b-(1-t)/t)/t^2
        s <- simpson_v2(f, 0, 1, n)
    } else if (a != -Inf & b == Inf) {
        f <- function(t) fun(a+(1-t)/t)/t^2
        s <- simpson_v2(f, 0, 1, n)
    } else {
        h <- (b-a)/n
        x <- seq(a, b, by=h)
        y <- fun(x)
        y[is.nan(y)]=0
        s <- y[1] + y[n+1] + 2*sum(y[seq(2,n,by=2)]) + 4 *sum(y[seq(3,n-1, by=2)])
        s <- s*h/3
    }
    return(s)
}

EX  <- function(x) x*(abs(x)* exp(-0.5*abs(x))+(1/(1+(x-40)^4)))
simpson_v2(EX, -Inf, Inf, n=100)

EX3 <- function(x) (x^3)*(abs(x)* exp(-0.5*abs(x))+(1/(1+(x-40)^4)))
simpson_v2(EX3, -Inf, Inf, n=100)

(ii) Importance Sampling My proposal density is Normal with mean=0 and standard deviation =4. The summary of the Importance sampling process I am applying is as follows

Suppose I can't sample from f(x) which is true as it has no well-known form and no built-in function is available in R to use for sampling. So, I propose another log cancave tail distribution N(0, 4) to take samples such that instead of estimating E(x) I estimate E(x*f(x)/N(0,1)). I use the following code for this which takes 100000 samples from N(0,4) 

X <- rnorm(1e5, sd=4)
Y <- (X)*(abs(X)*exp(-0.5*abs(X))+(1/(1+(x-40)^4)))/(dnorm(X, sd=4))
mean(Y)

Since this code needs random sampling from Normal distribution therefore, each time I got different answers but it is something around -0.1710694 which is almost similar to 0.3222642. I got from Simpsons method. But these results are very different E(X)= 88.85766 from integrate(). Note that integrate() use the Adaptive quadrature method. Is this method different from Simpsons and Importance sampling? What similarity in results I should expect while comparing these methods

Fatima Batool
  • 43
  • 7
  • 2
    This question is way too long. Also how do you propose to plot from -infinity to infinity? Just have an infinitely long x-axis? Do you want to do a transformation of some sort? – Dason Dec 31 '15 at 01:05
  • Yes, the question is too long but I have list every question separate so you can answer any specific part you wish to answer. But all are related to one problem. I am also confused with [-Inf, Inf] interval all the time. At first it makes no sense to plot a function over [-Inf, Inf]. But it's just a range and function can be plotted by longer gaps [-Inf, -1e20, -1e15, ...., +Inf]. If not should I just assume finite interval to approximate this integral. – Fatima Batool Dec 31 '15 at 01:20
  • No the target is not transformation. I just want to find E(x) and E(X^3) using various methods and try to compare them – Fatima Batool Dec 31 '15 at 01:22
  • I want to visualise the function over whole range. I think another idea is to identify the intervals where curve actually changes behaviour. And then find E(X) and E(X^3) for these intervals separately? – Fatima Batool Dec 31 '15 at 01:28
  • 1
    "I want to visualise the function over whole range". Once again how do you propose to incorporate an infinite range on a finite screen? And just because you broke it into different parts doesn't mean that's how questions are best handled here. You have a lot of non-programming related questions too. As it stands this "question" feels "too broad" and is in danger of being closed. – Dason Dec 31 '15 at 01:41
  • Ok let me try to narrow down the scope of the question – Fatima Batool Dec 31 '15 at 01:56
  • In general just ask a single easily defined problem. If it's more probability/statistics related then you should post those questions at http://stats.stackexchange.com/ instead of here. – Dason Dec 31 '15 at 01:58
  • Ohh I wasn't aware of this thank you for suggesting stats.stackexchange.com – Fatima Batool Dec 31 '15 at 02:03
  • Are you sure that `integrand` and `EX` is identical functions? `EX` has `exp(1/2x)` whereas `integrand` has `exp(-1/2x)` – Khashaa Dec 31 '15 at 02:58
  • 1
    `simpson_v2(integrand, 25, 50, n=1000)` is 88.843 – Khashaa Dec 31 '15 at 02:59
  • 1
    You failed to provide the probability density function for your random variable. Or you have to compute E(1) in your non-normalized version and divide by it to get the proper expectation values. – Lutz Lehmann Dec 31 '15 at 14:54
  • @LutzL Of course she provided the probability density function for your random variable. It is not normalized though, that's true, see my corrected answer – Severin Pappadeux Dec 31 '15 at 17:24
  • Can you rewrite your expression `f(x)∝ ...` with more parentheses so it's unambiguous? The way you've written it, it's equivalent to `(|x|exp(-x/2)) + (1/(1 + ((x-40)^4)))`. – Amit Kumar Gupta Dec 31 '15 at 23:03

1 Answers1

3

First, EX and EX3 definition is wrong, you miss minus under exponent

Well, here are some simplifications

UPDATE

Looks like your EX would be 40*\pi / \sqrt{2}

And EX3 is not infinity, I might be wrong here

UPDATE 2

Yep, EX3 is finite, should be a^2*EX + \pi*a*3/\sqrt{2}, where a is equal to 40

UPDATE 3

As noted, there is also a normalization required to get true values of EX and EX3

N = 8 + \pi/\sqrt{2}

Computed integrals to be divided by N to get proper moments.

Severin Pappadeux
  • 18,636
  • 3
  • 38
  • 64
  • upvote for analytic solution. `40*pi/sqrt(2)` is `88.85766` – Khashaa Dec 31 '15 at 05:30
  • @Khashaa thanks! Looks like `EX3` is good as well, `integrate` and analytical gives the same `142,438.8` answer – Severin Pappadeux Dec 31 '15 at 06:04
  • EX=5.2^(5/2)pi = 88.857658 Looks fine as Khashaa and Severin already mentioned. For EX3 from (http://www.integral-calculator.com/) I am getting graph which shows function have some behaviour from (-35, 0), (0, 25), (25, 55) and beyond these regions is almost zero. So basically there must be three bell curves in the region (-40, 60). When I am trying to plot this in R for (-35, 0), (0, 25), (25, 60) separately. I am getting one bell shape in each region but for (0, 55), (-35, 55) I am not getting two and three curves respectively for these regions. What I'am doing wrong? – Fatima Batool Dec 31 '15 at 20:42
  • the code is curve((x^3)*(abs(x)* exp(-0.5*abs(x))+(1/(1+(x-40)^4))), -35, 60) – Fatima Batool Dec 31 '15 at 20:44
  • Thanks @SeverinPappadeux for pointing minus sign mistake. – Fatima Batool Dec 31 '15 at 20:48
  • @SeverinPappadeux I am unable to handle Exp[-0.5abs(x)] in online mathematica? Which software you are using to calculate EX3= a^2*EX + \pi*a*3/\sqrt{2}? C?? – Fatima Batool Dec 31 '15 at 21:01
  • @FatimaBatool You don't have to calculate Exp() term at all. Why? Because it is same for `x` as well as for `-x`. So when you multiply it by odd power like `x` or `x^3` in the interval [-Inf...Inf] you'll get contribution from negative `x` term exactly cancelling contribution from positive `x` term. This is why I said that integrals from Exp() term for either `EX` or `EX3` are ALWAYS exactly ZERO, you could drop them from consideration even before going to online Mathematica tool. The only contribution from Exp() term is to the normalization. – Severin Pappadeux Dec 31 '15 at 21:17
  • thanks @SeverinPappadeux, what your suggestion about R graph problem? – Fatima Batool Dec 31 '15 at 21:48
  • @FatimaBatool It is hard to visualize what the problem is. What I would advice is to ask another pure plotting related question, no computation problem at al. Just show all of us function, plotting code, few plotting outputs and we will go from there. – Severin Pappadeux Dec 31 '15 at 22:14