0

I have recently been trying to convert a piece of Matlab code into Python code.

I have made most of the changes that I need to however, the issue I am having is the line where it says:

y(index(m)) = 1-x(index(m));

I get the error:

"Can't assign to function call"

However I am not sure how to restructure it in order to remove this error.

I have had a look around and people mention "get item" and "set item" however I have tried to use them, but I can't get them to work (probably because I can't figure out the structure)

Here is the full code:

import numpy

N = 100;
B = N+1;
M = 5e4;
burnin = M;
Niter = 20;
p = ones(B,Niter+1)/B;
hit = zeros(B,1);

for j in range(1,Niter):
    x = double(rand(1,N)>0.5);
    bin_x = 1+sum(x);
    index = ceil(N*rand(1,M+burnin));
    acceptval = rand(1,M+burnin);
    for m in range(1,M+burnin):
        y = x;
        y(index(m)) = 1-x(index(m));
        bin_y = 1+sum(y);

        alpha = min(1, p(bin_x,j)/p(bin_y,j) );
        if acceptval(m)<alpha:
            x = y; bin_x = bin_y;
        end

        if m > burnin: hit(bin_x) = hit(bin_x)+1; end
    end

    pnew = p[:,j];
    for b in range(1,B-1):
        if (hit(b+1)*hit(b) == 0):
            pnew(b+1) = pnew(b)*(p(b+1,j)/p(b,j));
        else:
            g(b,j) = hit(b+1)*hit(b) / (hit(b+1)+hit(b));
            g_hat(b) = g(b,j)/sum(g(b,arange(1,j)));
            pnew(b+1) = pnew(b)*(p(b+1,j)/p(b,j))+((hit(b+1)/hit(b))^g_hat(b));
        end
    end
    p[:,j+1] = pnew/sum(pnew);
    hit[:] = 0;
end

Thanks in advance

1 Answers1

3

The round brackets () indicate a function. For indexing you need [] square brackets - but that is only the first of many, many errors... I am currently going through line by line, but it's taking a while.

This code at least runs... you need to figure out whether the indexing is doing what you are expecting since Python arrays are indexed from zero, and Matlab arrays start at 1. I tried to fix that in a couple of places but didn't go through line by line - that's debugging.

Some key learnings:

  • There is no end statement... just stop indenting
  • When you import a library, you need to reference it (numpy.zeros, not zeros)
  • Lists are indexed from zero, not one
  • Indexing is done with [], not ()
  • Creating an array of random numbers is done with [random.random() for r in xrange(N)], not random(N).
  • ... and many other things you will find as you look through the code below.

Good luck!

import numpy
import random

N = int(100);
B = N+1;
M = 5e4;
burnin = M;
Niter = 20;
p = numpy.ones([B,Niter+1])/B;
hit = numpy.zeros([B,1]);
g = numpy.zeros([B, Niter]);
b_hat = numpy.zeros(B);

for j in range(1,Niter):
    x = [float(random.randint(0,1)>0.5) for r in xrange(N)];
    bin_x = 1+sum(x);
    index = [random.randint(0,N-1) for r in xrange(int(M+burnin))];
    #acceptval = rand(1,M+burnin);
    acceptval = [random.random() for r in xrange(int(M+burnin))];
    for m in range(1,int(M+burnin)):
        y = x;
        y[index[m]] = 1-x[index[m]];
        bin_y = 1+sum(y);

        alpha = min(1, p[bin_x,j]/p[bin_y,j] );
        if acceptval[m]<alpha:
            x = y; bin_x = bin_y;

        if m > burnin: 
            hit[bin_x] = hit[bin_x]+1;

    pnew = p[:,j];
    for b in range(1,B-1):
        if (hit[b+1]*hit[b] == 0):
            pnew[b+1] = pnew[b]*(p[b+1,j]/p[b,j]);
        else:
            g[b,j] = hit[b+1]*hit[b] / [hit[b+1]+hit[b]];
            g_hat[b] = g[b,j]/sum(g[b,numpy.arange(1,j)]);
            pnew[b+1] = pnew[b]*(p[b+1,j]/p[b,j])+((hit[b+1]/hit[b])^g_hat[b]);
    p[:,j+1] = pnew/sum(pnew);
    hit[:] = 0;
Floris
  • 45,857
  • 6
  • 70
  • 122
  • Thank you so much, you've been a great help, I'll go through it now and try and understand everything you've done. Thanks again –  Mar 13 '13 at 16:22