MATLAB conversion to python? Error in python due to shapes of arrays

Question

Question

I have a Brownian motion vectorized path given below that I am trying to replicate in python. My problem is that one of the functions, highlighted as the first function in yellow is not working. Instead it is giving me an error (operands could not be broadcast together with shapes (1000,499) (500,1000)) due to the dimensions of the two arrays being added.

How can I recreate the array U using python?

The MATLAB code is from An Algorithmic Introduction to Numerical Simulation of stochastic Differential Equations and is not my own code.

MATLAB code

Python attempt not working

import numpy as np
import matplotlib.pyplot as plt
from numpy import matlib as mb

# BPATH 3 Function along a Brownian path

T = 1
N = 500
dt = float(T/N)
M = 1000 # ntraj
t = np.arange(dt,T,dt)

## == Mean of 1000 paths == ##

dW = np.sqrt(dt)* np.random.normal(0.0, 1.0, (n,M)) # all general increments
W = np.cumsum(dW, axis=1) # cumsum by column

a = mb.repmat(t,M,1)
U = np.exp(a+0.5*W) # Error is here

Error message

     17 W = np.cumsum(dW, axis=1) # cumsum by column
     18 a = mb.repmat(t,M,1)
---> 19 U = np.exp(a+0.5*W)
     20 
     21 

ValueError: operands could not be broadcast together with shapes (1000,499) (500,1000)

Hi @wwii I have read the document but I am not understanding how this helps with the problem. Could you perhaps be more clear? — user4933
– user4933, Commented Nov 2, 2020 at 15:22
The document explains the rules for broadcasting. The array shapes need to be compatible for the operation being performed. If you understand the matlab code you wrote you should be able to craft a functional equivalent with numpy after you learn some numpy basics - that starts with spending some time with the numpy documentation. Unfortunately I am matlab challenged. — wwii
– wwii, Commented Nov 2, 2020 at 15:26
in the matlab code can you tell us the shape of each element? It will point the error relatively easy — jimakr
– jimakr, Commented Nov 2, 2020 at 15:47

hpaulj · Accepted Answer · 2020-11-02 19:09:49Z

2

You shouldn't need to import and use mb.repmat. That's used in MATLAB because it doesn't do broadcasting (or at least didn't until recently). To add to the (500,1000) shape W, t needs to be (500,), expended to (500,1). I'd suggest using linspace for generating t.

Test:

t = np.linspace(dt,T, N)   # (500,) shape
...
U = t[:,None] + 0.5*W      # add (500,1) with (500,1000)

np.arange doesn't include the end point so your t was only (499,) shaped.

===

The MATLAB code, reproduced in Octave

>> T=1;N=500;dt=T/N;t=[dt:dt:1];
>> M=1000;
>> x=repmat(t,[M 1]);
>> W=randn(M,N);
>> Umean=mean(x+W);

t is (1,500); x and W are (1000,500) size, and thus can add elementwise. Umean is (1,500) - the mean over the 1st dimension.

The numpy code:

In [1]: T = 1
   ...: N = 500
   ...: dt = float(T/N)
   ...: M = 1000 # ntraj
   ...: t = np.arange(dt,T,dt)
In [2]: t.shape
Out[2]: (499,)
In [3]: t = np.linspace(dt,T,N)
In [4]: t.shape
Out[4]: (500,)

In [6]: W=np.random.normal(0.0, 1.0, (N,M))
In [7]: W.shape
Out[7]: (500, 1000)

In [8]: (t[:,None]+W).shape
Out[8]: (500, 1000)

t[:,None] is (500,1) shape, which broadcasts to (500,1000) without using repmat.

To get the mean over the size 1000 dimension we have to use:

In [9]: np.mean((t[:,None]+W), axis=1).shape
Out[9]: (500,)

edited Nov 2, 2020 at 19:09

answered Nov 2, 2020 at 17:29

hpaulj

233k14 gold badges260 silver badges392 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

user4933 Over a year ago

Thanks @hpaulj . I see that the t is then (500, 1), should this not be (500, 1000) where 1000 represents M? When I run this the Umean that is being generated is non linear but the Umean is supposed to be a sloped line when you plot it rather than one with a lot of noise.

user4933 Over a year ago

The function elementwise would be like this: U[i,j] = np.exp(t+0.5*W[i])

hpaulj Over a year ago

I've added some Octave and numpy code that hopefully clarify the dimensions in both. Years ago I learned that getting the shapes right was 80% of the debugging battle in MATLAB. Don't guess; test and verify.

Collectives™ on Stack Overflow

MATLAB conversion to python? Error in python due to shapes of arrays

1 Answer 1

3 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

3 Comments

Your Answer

Sign up or log in

Post as a guest

Related