Matrix calculation (in R)  least squares
I have a matrix in R as per the following:
(72  90  60)
(60 83   )
(50  77 83 )
( 40  84 65)
Reproducible as follows:
B = matrix(c(72,60,50,NA,NA,83,NA,40,90,NA,77,NA,NA,NA,83,84,60,NA,NA,65),nrow=4,ncol=5)
It is a sparse matrix in the sense that not all elements are populated. What i want to do is find an adjustment matrix (or vector, that is added to each column) so that the average of the rows is equal to the average of columns. The answer is:
(13.04, 10.93, 14.68, 22.70, 5.90)
I computed this iteratively and is also rounded to 2dp so doesn't quite equate, but, for example, it should be the case that if you add the adjustments to the original matrix:
A = t(t(B)+c(13.04,10.93,14.68,21.70,5.9))
then the average of the row averages should equal the average of the column averages:
mean(colMeans(A, na.rm = TRUE, dims = 1))
mean(rowMeans(A, na.rm = TRUE, dims = 1))
However, I think this is a least squares problem and can be solved algebraicly but i don't know how to do it or how to set it out. The question is, is there a way to find the adjustment algebraicly using matrix calculation rather than iteratively. Can anyone help?
See also questions close to this topic

Mysql query works well at workbench but takes too long in r
I have a query to run in R which retrieves data from the database and performs operations on it. When I run it in mysql workbench, it works just fine but in r it takes way too long and may hang the entire system. I also tried to run it in command prompt but got the error:
Error: memory exhausted (limit reached?)
mysql query:
library(DBI) library(RMySQL) con < dbConnect(RMySQL::MySQL(), dbname ="mydb", host = "localhost", port = 3306, user = "root", password = "") pedigree < dbGetQuery (connection, "SELECT aa.name as person, mother as mom, father as dad FROM addweight LEFT JOIN aa ON addweight.name2 = aa.name2 or addweight.name = aa.name LEFT JOIN death ON addweight.name2 = death.name2 or addweight.name = death.name Where((death.dodeath > curdate() OR aa.name2 NOT IN (SELECT name2 FROM death) OR aa.name NOT IN (SELECT name FROM death) OR aa.name NOT IN (SELECT name FROM death)) AND (dob < curdate() AND domove < curdate()))")

Can I use Zeppelin as an alternative to Shiny?
I've read that Zeppelin can also do R visualizations using spark.r
My question is can I use it to do visualizations based on user inputs. These users would have not R/zeppelin technical experience.

R Write data in a file
I'm trying to save data in a file, but every time I hit the save button, it saves it but keeps deleting the data I already have there. What could be the problem?
saveData < function(data) { data < as.data.frame(t(data)) if (exists("responses")) { responses << rbind(responses, data) } else { responses << data } write.csv(responses, file = "read.csv", row.names = FALSE)

Why do transformations on a 3d box give random results on the screen?
I am having a problem where any transformations performed on a 3d model create random results. Most times the screen is blank without a gl error, sometimes it is like I am looking at the box from the inside, and on the few occasions code hacking has produced a 3d box the perspective it inverted where things farther away are larger. The orthographic projection works as intended so I believe it is an issue specific to the perspective matrix, but even after referencing perspective matrices several times I cannot find an issue with it. So that theory is in doubt.
Mat4.cpp (the constructer is in the header but it just set the matrix array to zero)
Mat4f Mat4f::Identity() { Mat4f result; result.m_Matrix[0 + 0 * 4] = 1; result.m_Matrix[1 + 1 * 4] = 1; result.m_Matrix[2 + 2 * 4] = 1; result.m_Matrix[3 + 3 * 4] = 1; return result; } Mat4f Mat4f::translate(const Vec3f& vector) { Mat4f result = Identity(); result.m_Matrix[0 + 3 * 4] = vector.m_X; result.m_Matrix[1 + 3 * 4] = vector.m_Y; result.m_Matrix[2 + 3 * 4] = vector.m_Z; return result; } Mat4f Mat4f::rotate(const Vec3f& axis, const float& angle) { Mat4f result = Identity(); float radians = toRadians(angle); float cosResult = std::cos(radians); float sinResult = std::sin(radians); float cosInv = 1.0f  cosResult; float x = axis.m_X; float y = axis.m_Y; float z = axis.m_Z; //first column result.m_Matrix[0 + 0 * 4] = x * cosInv + cosResult; result.m_Matrix[1 + 0 * 4] = sinResult * z; result.m_Matrix[2 + 0 * 4] = sinResult * y; //second column result.m_Matrix[0 + 1 * 4] = sinResult * z; result.m_Matrix[1 + 1 * 4] = y * cosInv + cosResult; result.m_Matrix[2 + 1 * 4] = sinResult * x; //third column result.m_Matrix[0 + 2 * 4] = sinResult * y; result.m_Matrix[1 + 2 * 4] = sinResult * x; result.m_Matrix[2 + 2 * 4] = z * cosInv + cosResult; return result; } Mat4f& Mat4f::multiply(const Mat4f& matrix) { Mat4f result; for (int column = 0; column < MAT4F_DIMINSION; column++) { for (int row = 0; row < MAT4F_DIMINSION; row++) { for (int element = 0; element < MAT4F_DIMINSION; element++) { result.m_Matrix[column + row * MAT4F_DIMINSION] += m_Matrix[column + element * MAT4F_DIMINSION] * matrix.m_Matrix[element + row * MAT4F_DIMINSION]; } } } m_Matrix.swap(result.m_Matrix); return *this; } Mat4f& Mat4f::operator*=(const Mat4f& matrix) { return multiply(matrix); } Mat4f operator*(Mat4f& left, const Mat4f& right) { return left.multiply(right); } Mat4f Mat4f::Perspective(float aspectRatio, float fov, float far, float near) { Mat4f result; float normalizedFOV = toRadians(fov); float tanResult = std::tan(normalizedFOV / 2); result.m_Matrix[0 + 0 * 4] = (1 / (tanResult)) /aspectRatio; result.m_Matrix[1 + 1 * 4] = 1 / tanResult; result.m_Matrix[2 + 2 * 4] = 1 / (far + near) / (near  far); result.m_Matrix[3 + 2 * 4] = 1; result.m_Matrix[2 + 3 * 4] = 1 / (2 * far * near) / (near  far); return result; }
main.cpp (the vertices in this file are set up so that the first numbers on each line are positions and the next four is color data)
using namespace Engine; int main() { Engine::Window* window = Engine::Window::createWindow(1280, 720, "test"); Graphics::Shader shader("res/Shaders/BasicVert.sdr", "res/Shaders/BasicFrag.sdr"); shader.bind(); float vertices[] = { //front 1.0, 1.0, 1.0, 1.0f, 1.0f, 1.0f, 1.0f, 1.0, 1.0, 1.0, 0.0f, 1.0f, 1.0f, 1.0f, 1.0, 1.0, 1.0, 1.0f, 0.0f, 1.0f, 1.0f, 1.0, 1.0, 1.0, 1.0f, 1.0f, 0.0f, 1.0f, // back 1.0, 1.0, 1.0, 0.0f, 0.0f, 1.0f, 1.0f, 1.0, 1.0, 1.0, 1.0f, 0.0f, 0.0f, 1.0f, 1.0, 1.0, 1.0, 0.0f, 1.0f, 0.0f, 1.0f, 1.0, 1.0, 1.0, 0.0f, 0.0f, 0.3f, 1.0f, }; unsigned int indecies[] = { 0, 1, 2, 2, 3, 0, // top 1, 5, 6, 6, 2, 1, // back 7, 6, 5, 5, 4, 7, // bottom 4, 0, 3, 3, 7, 4, // left 4, 5, 1, 1, 0, 4, // right 3, 2, 6, 6, 7, 3, }; Math::Mat4f projection = Math::Mat4f::Perspective(((float)16/9), 45, 10, .1); Math::Mat4f rotation = Math::Mat4f::rotate({ 0, 1, 0}, 45); Math::Mat4f transform = Math::Mat4f::Identity(); transform *= rotation; transform *= Math::Mat4f::translate({0, 0, 1}); Graphics::VertexBufferLayout layout; layout.pushFloats(3); layout.pushFloats(4); Graphics::VertexArray va; Engine::Graphics::Buffer buffer(vertices, sizeof(vertices)); buffer.bind(); Engine::Graphics::IndexBuffer indexbuffer(indecies, 36); indexbuffer.bind(); va.addBuffer(buffer, layout); glEnable(GL_DEPTH_TEST); while (!glfwWindowShouldClose(window>m_Windowptr)) { glClear(GL_COLOR_BUFFER_BIT  GL_DEPTH_BUFFER_BIT); shader.setUniformMat4("m_Projection", projection); shader.setUniformMat4("m_Transform", transform); std::cout << glGetError() << std::endl; glDrawElements(GL_TRIANGLES, 36, GL_UNSIGNED_INT, 0); glfwSwapBuffers(window>m_Windowptr); glfwPollEvents(); } }
Vertex Shader
#version 330 core layout(location = 0)in vec4 position; layout(location = 1)in vec4 in_color; uniform mat4 m_Projection; uniform mat4 m_Transform; out vec4 vs_color; void main(){ gl_Position = m_Projection * m_Transform * position; vs_color = in_color; }
Let me know if you need any other file.

Neural Network loss function and weight matrix doesn't match
I'm attempting to learn about neural networks using matrices and in doing so I decided to set myself a challenge, train a simple neural network to output 0.5*sigma(x).
My guess is this would be sufficient as W needs to be 0.5.
But I've run into a problem regarding the loss function and the weights.
import sys import math import numpy as np import random @np.vectorize def dlossbydw(X,W,Y): t = np.dot(X,W) Yhat = sig(t) n = YYhat l = 2*n * 1 * sig(t)*(1sig(t))*X return l @np.vectorize def sig(Z): return (1.0/ (1.0 + math.exp(Z))) @np.vectorize def toMatch(Z): return 0.5*sig(Z) def main(args): #This is the matrix that should be output Y = np.array([[toMatch(x)] for x in range(5,5)]) #This is the input matix X = np.array([[x] for x in range(5,5)]) random.seed(5) r = random.random() #And this is the weight matrix W = np.array([r]) rate = 1e1 for i in range(1000): print("dlossbydw: " + str(dlossbydw(X,W,Y))) #??? print("expected out:" + str(Y)) print("post training:" + str(sig(np.dot(X,W)))) if(__name__=="__main__"): main(sys.argv[1:])
The dimensions of the weight matrix W is 1x1, (there is only one weight connecting the input neuron and the output) however when I do the math to calculate the loss with respect to w, it gives me a 10x1 matrix, my guess is it gives me the loss for each X inputted for the input matrix, but what do I do with this? Each value in the loss is different from the last, which is strange considering they all require the same shift in W. (towards 0.5)
I must be misunderstanding something, or confusing myself over a simple mistake.
Can someone please clarify what I'm doing wrong here?
What should I be doing at the place marked?
Is this correct?
Thank you.

Create Correlation Matrix
I'd like to use the distance correlation matrix form the energy package to create all pairwise correlations, like
cor
would do:cor(iris[,1:4]) Sepal.Length Sepal.Width Petal.Length Petal.Width Sepal.Length 1.0000000 0.1175698 0.8717538 0.8179411 Sepal.Width 0.1175698 1.0000000 0.4284401 0.3661259 Petal.Length 0.8717538 0.4284401 1.0000000 0.9628654 Petal.Width 0.8179411 0.3661259 0.9628654 1.0000000
How can I do the same with
energy::dcor()
?> energy::dcor(iris[,1], iris[,2], 1.0) [1] 0.3105326 > energy::dcor(iris[,1], iris[,3], 1.0) [1] 0.8585197 > energy::dcor(iris[,1], iris[,4], 1.0) [1] 0.8266021 > ##etc...

About adding regularization to to the cost function of neural network for regression in tensorflow
I am trying to build a neural network for linear rgression. I want to add regularization part to the cost function but the cost does not change after each iteration. The code is as follows:
X = tf.placeholder(tf.float32,[n_x, None], name = "x") Y = tf.placeholder(tf.float32,[n_y, None], name = "y") W1 = tf.get_variable("W1", [25,11], initializer = tf.contrib.layers.xavier_initializer()) b1 = tf.get_variable("b1", [25,1], initializer = tf.zeros_initializer()) W2 = tf.get_variable("W2", [25,25], initializer = tf.contrib.layers.xavier_initializer()) b2 = tf.get_variable("b2", [25,1], initializer = tf.zeros_initializer()) W3 = tf.get_variable("W3", [25,25], initializer = tf.contrib.layers.xavier_initializer()) b3 = tf.get_variable("b3", [25,1], initializer = tf.zeros_initializer()) W4 = tf.get_variable("W4", [25,25], initializer = tf.contrib.layers.xavier_initializer()) b4 = tf.get_variable("b4", [25,1], initializer = tf.zeros_initializer()) W5 = tf.get_variable("W5", [12,25], initializer = tf.contrib.layers.xavier_initializer()) b5 = tf.get_variable("b5", [12,1], initializer = tf.zeros_initializer()) W6 = tf.get_variable("W6", [1,12], initializer = tf.contrib.layers.xavier_initializer()) b6 = tf.get_variable("b6", [1,1], initializer = tf.zeros_initializer()) Z1 = tf.add(tf.matmul(W1,X),b1) A1 = tf.nn.relu(Z1) Z2 = tf.add(tf.matmul(W2,A1),b2) A2 = tf.nn.relu(Z2) Z3 = tf.add(tf.matmul(W3,A2),b3) A3 = tf.nn.relu(Z3) Z4 = tf.add(tf.matmul(W4,A3),b4) A4 = tf.nn.relu(Z4) Z5 = tf.add(tf.matmul(W5,A4),b5) A5 = tf.nn.relu(Z5) Z6 = tf.add(tf.matmul(W6,A5),b6) A6 = tf.nn.tanh(Z6) regulaizers = tf.nn.l2_loss(W1) + tf.nn.l2_loss(W2) + tf.nn.l2_loss(W3) + tf.nn.l2_loss(W4) + tf.nn.l2_loss(W5) + tf.nn.l2_loss(W6) beta=0.01 cost = (1/(2*m))*tf.reduce_sum(tf.pow(A6Y, 2) ) + beta * regulaizers optimizer = tf.train.AdamOptimizer(learning_rate = learning_rate).minimize(cost) epoch=1500 init = tf.global_variables_initializer() with tf.Session() as sess: sess.run(init) for epoch in range(num_epochs): epoch_cost = 0. num_minibatches = int(m / minibatch_size) minibatches = random_mini_batches(X_train, Y_train, minibatch_size, seed) for minibatch in minibatches: (minibatch_X, minibatch_Y) = minibatch _ , minibatch_cost = sess.run([optimizer, cost], feed_dict={X: minibatch_X, Y: minibatch_Y}) epoch_cost += minibatch_cost / num_minibatches
After runnining the script through epochs, the output is
Cost after epoch 0: 4.524896 Cost after epoch 100: 0.000041 Cost after epoch 200: 0.000041 Cost after epoch 300: 0.000041 Cost after epoch 400: 0.000041 Cost after epoch 500: 0.000041 Cost after epoch 600: 0.000041 Cost after epoch 700: 0.000041 Cost after epoch 800: 0.000041 Cost after epoch 900: 0.000041 Cost after epoch 1000: 0.000041 Cost after epoch 1100: 0.000041 Cost after epoch 1200: 0.000041 Cost after epoch 1300: 0.000041 Cost after epoch 1400: 0.000041 Cost after epoch 1500: 0.000041 Cost after epoch 1600: 0.000041 Cost after epoch 1700: 0.000041 Cost after epoch 1800: 0.000041 Cost after epoch 1900: 0.000041
After initialising the parameter and feed the session, As you can see the cost does not change. I was wondering if I can get some help and if the cost function is correct.

Controlling for i.X and X produce the opposite result in Stata
I did the following
reg Y i.X, nocons coefplot, vertical keep(*.X)
and this produced a graph of very steeply decreasing pattern at almost every point.
However,
reg Y X, nocons
produced a statistically very significant positive coefficient on X!
In what cases would it be possible?

Why `sklearn` and `statsmodels` implementation of OLS regression give different R^2?
Accidentally I have noticed, that OLS models implemented by
sklearn
andstatsmodels
yield different values of R^2 when not fitting intercept. Otherwise they seems to work fine. The following code yields:import numpy as np import sklearn import statsmodels import sklearn.linear_model as sl import statsmodels.api as sm N=1000 X = np.random.normal(loc=1, size=(N, 1)) Y = 2 * X.flatten() + 4 + np.random.normal(size=N) sklernIntercept=sl.LinearRegression(fit_intercept=True).fit(X, Y) sklernNoIntercept=sl.LinearRegression(fit_intercept=False).fit(X, Y) statsmodelsIntercept = sm.OLS(Y, sm.add_constant(X)) statsmodelsNoIntercept = sm.OLS(Y, X) print(sklernIntercept.score(X, Y), statsmodelsIntercept.fit().rsquared) print(sklernNoIntercept.score(X, Y), statsmodelsNoIntercept.fit().rsquared) print(sklearn.__version__, statsmodels.__version__)
prints:
0.808887097312 0.808887097312 0.724295712882 0.786022399826 0.18.1 0.8.0
Where the difference comes from?
The question differs from Different Linear Regression Coefficients with statsmodels and sklearn as there
sklearn.linear_model.LinearModel
(with intercept) was fit for X prepared as forstatsmodels.api.OLS
. 
scipy.optimize.least_squares  limit number of jacobian evaluations
I am trying to use scipy.optimize.least_squares(fun= my_fun, jac=my_jac, max_nfev= 1000) with two callable functions: my_fun and my_jac
both fuctions: my_fun and my_jac, use an external software to evaluate their value, this task is much time consuming, therefore I prefer to control the number of evaluations for both
the trf method uses the my_fun function for evaluating if trust region is adequate and the my_jac function for determine the step direction
There is an input parameter max_nfev. does this parameter count only for the fun evaluations? does it consider also the jac evaluations?
moreover, in matlab there are two parameters for the lsqnonlin function, MaxIterations and MaxFunctionEvaluations. does it exist in scipy.optimize.least_squares?
Thanks Alon

R Least Squares Optimization Best Methods
I have a set of data, the production of oil from a well. And an equation to predict that forecast.
The equation requires 5 input variables which are real numbers with upper and lower bounds, 1 input variable which must be an integer and 1 input variable which can be 1 of 2 string variables.
What is the best optimization routine (I am using nlminb) to use in R to determine the best set of inputs and how to I handle the integer and string inputs when optimizing?
I have attached a plot of what I am trying to accomplish as well as my sample code.
Raw Data & Best Fit
sse < function(q, forecast) { sum(na.omit((q  forecast)) ^ 2) } bounds.low < function(Qi.low, b.low, Db.low, Df.low, Tper, Qo.low, Td.low) { if (Tper == "D"  Tper == "d") { Tper < 365 } else if (Tper == "M"  Tper == "m") { Tper < 12 Qi.low < Qi.low * (365 / 12) Qo.low < Qo.low * (365 / 12) } else if (Tper == "A"  Tper == "a") { Tper < 1 Qi.low < Qi.low * 365 Qo.low < Qo.low * 365 } c(Qi.low, b.low, Db.low, Df.low, Qo.low, Td.low) } bounds.high < function(Qi.high, b.high, Db.high, Df.high, Tper, Qo.high, Td.high) { if (Tper == "D"  Tper == "d") { Tper < 365 } else if (Tper == "M"  Tper == "m") { Tper < 12 Qi.high < Qi.high * (365 / 12) Qo.high < Qo.high * (365 / 12) } else if (Tper == "A"  Tper == "a") { Tper < 1 Qi.high < Qi.high * 365 Qo.high < Qo.high * 365 } c(Qi.high, b.high, Db.high, Df.high, Qo.high, Td.high) } best.modhypbu.vol < function(vol, t, lower, upper, Qf, Tmax, Tper, build.method) { # if (length(q) != length(t)  length(q) <= 2) # stop("Invalid lengths for q, t vectors.") #convert Units: if (Tper == "D"  Tper == "d") { Tper < 365 Tmax < Tmax * Tper } else if (Tper == "M"  Tper == "m") { Tper < 12 Qf < Qf * (365 / 12) Tmax < Tmax * Tper } else if (Tper == "A"  Tper == "a") { Tper < 1 Qf < Qf * 365 Tmax  Tmax * Tper } res < nlminb( c(# initial guesses max(na.omit(vol)), # qi 1.2, # b = rightish for a lot of wells currently coming on 0.75, # Di = decline from first to second point 0.06, # Df = about 9% effective vol[1], # Qo match(max(na.omit(vol)),vol)), # Td # cost function function (guess) sse( vol, modhypbu.vol( guess[1], guess[2], guess[3], guess[4], Qf, Tmax, Tper, t, guess[5], guess[6], build.method ) ), # bounds lower = lower, upper = upper ) ssr < sum(na.omit((vol  mean(na.omit( vol )))) ^ 2) sst < res$objective + ssr R2 < ssr / sst list( Qi = res$par[1] / (365 / Tper), b = res$par[2], Db = res$par[3], Df = res$par[4], Qo = res$par[5], Td = res$par[6], sse = res$objective, R2 = R2 ) } vol < DATA$OIL[1:20] #vol < oil #t < DATA$MONTH[2:20] t < (1:length(vol)) Qi.low < 1 b.low < 0 Db.low < 0.1 Df.low < 0.06 Qi.high < 2000 b.high < 2 Db.high < 1 Df.high < 0.06 Qo.low < 1 Qo.high < 50 Td.low < 1 Td.high < 1 build.method < "Exp" Qf < 1 Tmax < 50 Tper < "M" test < best.modhypbu.vol( vol, t, bounds.low(Qi.low, b.low, Db.low, Df.low, Tper, Qo.low, Td.low), bounds.high(Qi.high, b.high, Db.high, Df.high, Tper, Qo.high, Td.high), Qf, Tmax, Tper, build.method ) forecast < modhypbu.vol(test$Qi, test$b, test$Db, test$Df, Qf, Tmax, Tper, t, test$Qo, test$Td, build.method) plot(t, vol, type = "l", col = "red") lines(t, forecast, col = "green") test

Storing coefficients from a Regression in Stata
I am trying to store the coefficients from a simulated regression in a variable b1 and b2 in the code below, but I'm not quite sure how to go about this. I've tried using
return scalar b1 = _b[x1]
andreturn scalar b2 = _b[x2]
, from therclass()
function, but that didn't work. Then I tried usingscalar b1 = e(x1)
andscalar b2 = e(x2)
, from theeclass()
function and also wasn't successful.The goal is to use these stored coefficients to estimate some value (say rhat) and test the standard error of rhat.
Here's my code below:
program montecarlo2, eclass clear version 11 drop _all set obs 20 gen x1 = rchi2(4)  4 gen x2 = (runiform(1,2) + 3.5)^2 gen u = 0.3*rnormal(0,25) + 0.7*rnormal(0,5) gen y = 1.3*x1 + 0.7*x2 + 0.5*u * OLS Model regress y x1 x2 scalar b1 = e(x1) scalar b2 = e(x2) end
I want to do something like,
rhat = b1 + b2, and then test the standard error of rhat.