Computational Neuroscience Homework 1 Solution

Instructions:




Prepare a report (including your answers/plots) to be uploaded on Moodle.



The report should be typeset (no handwriting allowed except for lengthy derivations, which may be scanned and embedded into the report).



Show all the steps of your work clearly.



Unclear presentation of results will be penalized heavily.



No partial credits to unjusti ed answers.



Use Matlab or Python for computations.



Return all Matlab/Python code that you wrote in a single le.



Code should be commented, code for di erent HW questions should be clearly sepa-rated.



The code le should NOT return an error during runtime.



If the code returns an error at any point, the remaining part of your code will not be evaluated (i.e., 0 points).





















Question Points Your Score







Q1 25







Q2 25







Q3 25







Q4 25







TOTAL 100




Question 1. [25 points] Assume that a neural population computes weighted linear combinations of its input x, characterized by a system of equations Ax = b. Here A is the transfer function and b is the output vector.






0
1
0


1
2
1
x1




0
1
1




x3












B
x


C








A single output measurement is recorded, given by


3
3
0
9


=


9


.
@
2
1


1
5
0
x2
1
@
4
A












A B


4
C




















@




A










Answer the questions below. Derive the results by hand rst and then con rm them on the computer.




Find all solutions xn to Ax = 0.



Find a particular solution xp to Ax = b.



Find all solutions to Ax = b.



Find the pseudo-inverse of A.



Find the sparsest solution to Ax = b (the solution with the least number of non-zero entries).



Find the least-norm solution to Ax = b (the solution with the minimum Euclidean norm).
 
Question 2. [25 points]



We record the aggragated responses of a population of 100,000 neurons (rough number of neurons contained in a volumetric-pixel of fMRI acquisitions). 30 separate recordings of the aggragate response X are given in Xdata.mat. Answer the questions below. Derive the results on the computer.







Find the sample mean, and the sample median.



Find the sample standard deviation, and the sample inter-quartile range.



Compute the histogram of the data for the value range 70 X 130. Plot the histograms for n = 3, 6, 12 bins.



Determine whether X is normally distributed by running the following code on samples of X. (This is called the Normal Quantile Plot)



y = sort(x);




n = length(x);




f = ((1:n)-3/8)/(n+1/4);




q = 4.91*(f.^0.14 - (1-f).^0.14);




figure(1);clf;plot(q,y,'*-');grid;




Using bootstrap resampling with N = 1000 iterations, determine the sample mean and its standard error. Find the 95% con dence interval of the mean. Plot a histogram of the bootstrap distribution for n = 50 bins.



Using bootstrap resampling with N = 1000 iterations, determine the sample standard deviation and its standard error. Find the 95% con dence interval of the standard deviation. Plot a histogram of the bootstrap distribution for n = 50 bins.



Repeat parts e and f using jackknife resampling. Note that for jackknife, you resample by leaving one-sample-out at each iteration, so N = 30. Plot histograms for n = 5 bins.



Question 3. [25 points]




`Reverse inference' is a common, albeit poorly exercised method in neuroscience. It refers to the practice of inferring that a cognitive process is engaged on the basis of activation in some brain area. For example, if Broca's area was found to be activated in some task, researchers might infer that the subjects were using language. After a comprehensive search of the literature, we nd that Broca's area was reported to be activated in 103 out of 869 fMRI tasks involving engagement of language, but this area was also active in 199 out of 2353 tasks not involving language.







Assume that the conditional probability of activation given language and activation given no language, each follow a Bernoulli distribution (i.e., active with some probability p, or



not with probability 1 p). Compute the likelihoods of observed frequencies of activation in literature, as functions of the possible values of their respective Bernoulli probability pa-rameters p = xl and p = xnl. Compute these functions at the values x =[0:.001:1] and plot them as separate bar charts.




Find the values of xl and xnl that maximize their respective discretized likelihood func-tions.



Using the likelihood functions computed for discrete x, compute and plot the discrete pos-



terior distributions P (Xjdata) and the associated cumulative distributions P (X xjdata) for both processes (language and no language cases). To do this, assume a uniform prior




P (x) / 1 and note that it will be necessary to compute (rather than ignore) the normalizing constant for Bayes' rule. Use the cumulative distributions to compute (discrete approxima-tions to) upper and lower 95% con dence bounds on each proportion (xl;nl).




Consider the joint posterior distribution P (Xl; Xnljdata) over xl and xnl, the Bernoulli probability parameters for the language and non-language contrasts. Given that these two frequencies are independent, the (discrete) joint distribution is given by the outer product of



the two marginals. Plot it (with imagesc). Compute (by summing the appropriate entries in the joint distribution) the posterior probabilities that P (Xl Xnljdata) and conversely that P (Xl Xnljdata).




Using the estimates from part b as the relevant conditional probabilities, and assuming the prior that a contrast engages language, P (language) = 0:5, compute the probability



P (languagejactivation) that observing activation in this area implies engagement of lan-guage processes. Is the critique on `reverse inference' correct? How con dent should you be in implicating language if you observe activity in Broca's area?Question 4. [25 points]




Write a function samples = ndRandn(mean, cov, num) that generates a set of samples drawn from a multidimensional Gaussian distribution with the speci ed mean (an N-vector) and covariance (an NxN matrix). The parameter num should be optional (defaulting to 1) and should specify the number of samples to return. The returned value should be a matrix with num rows each containing a sample of N elements.







Your function should use the MATLAB function randn to generate samples from an N-



dimensional Gaussian with identity covariance matrix, and then modify these appropriately. (Hint: Recall that the sample covariance is 1=N XT X. If 1=N XT X is I then (XY )T (XY )

has sample covariance 1=N Y T Y . Therefore you need to multiply samples X by a matrix

such that Y T Y = cov. You can nd such a matrix using the SVD of cov).



Test your function by plotting 1000 samples of a 2-dimensional Gaussian (choose an arbitrary nonzero mean and nonzero covariance). Measure the sample mean and covariance of your data points, comparing to the values that you requested when calling the function. Plot an ellipse on top of the scatterplot that traces out points that are two standard devi-ations away from the mean, according to the covariance matrix. Does this ellipse capture the shape of the data?