Sumation inside a for loop

General Usage
1

i have tried to do a simple code for MCMC as you see below:

traget(x)= (1/sqrt(2*pi*2^2))*exp((-(x-10)^2)/(2*2^2))
proposal(x)=Normal(x,1)
initial=-40
result=[]
Iteration=10000
chain=zeros(Iteration)
accept_rate= zeros(Iteration)

for i in 1:Iteration
    accept=0
    current= initial
    proposed= rand(proposal(current)) 
    C= traget(proposed)/traget(current)
    ratio= min(1, C)
    uniform=rand()
    if uniform< ratio
        push!(result, proposed)
        current= proposed
        
        accept+=1
    else
        current=current
        
    end
    chain[i]= current
    accept_rate[i]= accept/i
end

my problem is when i want to know how many values have passed the test which i can know by length(result) or just accept but when i type accept it shows me only either 0 or 1 and i want to show the number of the proposed values that have passed the MCMC test

2

You reset accept to zero at every step of your iteration, therefore it can only ever be 1 (if the last iteration produced a proposal that was accepted) or 0 (if it didn’t).

You prpbably meant to set accept = 0 outside you iteration loop.

1 Like
3

If the above code runs in the REPL and if you put accept = 0 outside the loop, don’t forget to change the counting line to

global accept+=1

Looking at your code, I guess, that you put accept=0 into the loop, because you run into some error messages you didn’t fully understand.
If this is the case, you may have some insights reading this:
https://docs.julialang.org/en/v1/manual/variables-and-scoping/

3 Likes
4

In addition to the advice above: it is much easier to work with this kind of code in functions instead of “scripts”.

Incidentally, you almost surely want to work with log pdfs.

2 Likes
5

I am afraid this needs some further explanation :wink:

6

Log probability density functions.

log_target(x) = (-(x-10)^2)/(2*2^2)

...

log_C = log_target(proposed) - target(current)

etc.

2 Likes
7

= pdf

I was stuck in mind with Portable Document Format :smiley:

8

just to make sure you mean log_target(proposed)-log_target(current)right ? not log_target(proposed) - target(current)

1 Like
9

Sure, that’s a typo.