text:
# Difference between revisions of "sparse PCA"

## Introduction

In PCA, Given $n$ observations on $d$ variables (or in other words $n$ $d$-dimensional data points), our goal is to find directions in the space of the data set that correspond to the directions with biggest variance in the input data. In practice each of the $d$ variables has its own special meaning and it may be desirable to come up with some directions, as principal components, each of which is a combination of just a few of these variables. This makes the directions more interpretable and meaningful. But this is not something that usually happens as the original result of PCA method. Each of resulting directions from PCA in most cases is a linear combination of all variable with no zero coefficients.

To address the above concerns we add a sparsity constraint to the PCA problem, which makes the PCA problem much harder to solve. That's because we have just added a combinatorial constraint to optimization problem. This paper is showing us how to find directions in the data space with maximum variance that have a limited number of non-zero elements. In other words, this helps us to perform feature selection, by selecting a subset of features in each direction.

## Contribution

In this paper, a direct approach (called DSPCA) that improves the sparsity of the principle components is presented. This is done in 2 stages. First, incorporating a sparsity criterion in the PCA formulation. Second, forming a convex relation of the problem that is a semidefinite program. For small problems, semidifinite programs can be solved via general purpose interior-point methods. However, these methods can not be used for high dimensional problems. In this case, a saddle point problem can express our particular problem. For this kind of problems, smoothing argument algorithms combined with an optimal first-order smooth minimization algorithm offer a significant reduction in computational time and therefore can be used instead of generic interior point SDP solvers.

## Notation

The following notations are used in this note.

$S^n \,$ is the set of symmetric matrices of size $n \,$.
$\textbf{1} \,$ is a column vector of ones.
$\textbf{Card}(x) \,$ denotes the cardinality (number of non-zero elements) of a vector $x \,$
$\textbf{Card}(X) \,$ denotes the cardinality (number of non-zero elements) of a matrix $X \,$
For $X \in S^n \,$, $X \succeq 0 \,$ means $X \,$ is positive semi-definite.
$|X| \,$ is the matrix whose elements are the absolute values of the elements of $X \,$

## Problem Formulation

Given the covariance matrix $A$, the problem can be written as:

 $\begin{array}{ll} \textrm{maximize}& x^TAx\\ \textrm{subject\ to}& ||x||_2=1\\ &\textbf{Card}(x)\leq k \end{array}$ (1)

The cardinality constraint makes this problem hard (NP-hard) and we are looking for a convex and efficient relaxation.

Defining $X=x^Tx$, the above formula can be rewritten as

 $\begin{array}{ll} \textrm{maximize}& \textbf{Tr}(AX)\\ \textrm{subject\ to}&\textbf{Tr}(X)=1\\ &\textbf{Card}(X)\leq k^2\\ &X\succeq 0, \textbf{Rank}(X)=1\\ \end{array}$ (2)

The conditions $X\succeq 0$ and $\textbf{Rank}(X)=1$ in formula 2 guarantees that $X$ can be written as $x^Tx$, for some $x$. But this formulation should be relaxed before it can be solved efficiently, because the constraintS $\textbf{Card}(X)\leq k^2$ and $\textbf{Rank}(X)=1$ are not convex. So we replace the cardinality constraint with a weaker one: $\textbf{1}^T|X|\textbf 1\leq k$. We also drop the rank constraint. So we get:

 $\begin{array}{ll} \textrm{maximize}& \textbf{Tr}(AX)\\ \textrm{subject\ to}&\textbf{Tr}(X)=1\\ &\textbf{1}^T|X|\textbf 1\leq k\\ &X\succeq 0\\ \end{array}$

The above semidefinite relaxation can even be generalised to a non square matrix $A \in R^{mxn}$ as follows:

 $\begin{array}{ll} \textrm{maximize}& \textbf{Tr}(AX^{12})\\ \textrm{subject\ to}&\textbf{Tr}(X^{ii})=1\\ &\textbf{1}^T|X^{ii}|\textbf 1\leq k_i, i=1,2\\ &\textbf{1}^T|X^{12}|\textbf 1\leq \sqrt{k_1k_2}\\ &X\succeq 0\\ \end{array}$

in the variable $X \in S^{m+n}$ with blocks $X^{ij}$ for i,j=1,2.

We then change the modified cardinality constraint to a penalty term in the goal function with some positive factor $\rho$. So we get a semidefinite form of the problem:

 $\begin{array}{ll} \textrm{maximize}& \textbf{Tr}(AX)-\rho\textbf{1}^T|X|\textbf 1\\ \textrm{subject\ to}&\textbf{Tr}(X)=1\\ &X\succeq 0\\ \end{array}$ (3)

where, $\rho$ controls the penalty magnitude.

The goal function can be rewritten as $\textbf{Tr}(AX)-\rho\textbf{1}^T|X|\textbf 1=\min_{|U_{ij}|\leq\rho}\textbf{Tr}((A+U)X)$. So the problem (3) is equivalent to:

 $\begin{array}{ll} \textrm{maximize}& \min_{|U_{ij}|\leq\rho}\textbf{Tr}(X(A+U))\\ \textrm{subject\ to}&\textbf{Tr}(X)=1\\ &X\succeq 0\\ \end{array}$ (4)

or equivalently, due to convexity:

 $\begin{array}{lll} \textrm{minimize}& \lambda^{\max}(A+U)\\ \textrm{subject\ to}&|U_{ij}|\leq\rho,&i,j=1\cdots,n\\ \end{array}$ (5)

where $\lambda^{\max}(M)$ is the largest eigenvalue of the matrix $M$.

The problem as described in formulation (5) can be seen as computing a robust version of maximum eigenvalue: it is the least possible value of maximum eigenvalue, given that each element can be changed by at most noise value $\rho$. Also, it corresponds to a worst-case maximum eigenvalue computation with a bounded noise of intensity $\rho$ in each component on the matrix coefficients.

The KKT conditions for optimization problems (3) and (5) are given by:

 $\left\{ \begin{array}{rl} &(A+U)X=\lambda^{\max}(A+U)X\\ &U\circ X=\rho |X| \\ &\text{Tr}(X)=1,\,\,\,X\succeq 0 \\ &|U_{i,j}|\leq \rho,\,\,\, i,j=1,\cdots ,n \end{array} \right.$

If the $\lambda^{\max}$ in the first equation is simple (meaning it is of multiplicity 1) and $\rho$ is sufficiently small, from the first equation it follows that $\textbf{Rank}(X)=1$. In fact, the form of this equation implies that all columns of $\,X$ are eigenvectors of matrix $A+U$ corresponding to its maximum eigenvalue. So the rank one constraint is automatically satisfied in this special case.

## The Algorithm

### The Main Loop

The algorithm should iteratively create the semidefinite program (4) and solve it to obtain the next most important sparse principle component. At each iteration, we should first obtain the solution $x$ of the corresponding problem of form (1), if $X$ is the optimal solution of the optimization problem. That will be straightforward if $X$ is of rank 1, but since we have dropped the rank constraint, this may not be true and in those cases we need to obtain the dominant eigenvalue of $X$ by the methods that are known in the literature; for example we can use the power method which efficiently provides us with the largest eigenvectors of a matrix. Note, however, that in this case the resulting vectors are not guaranteed to be as sparse as the matrix itself. After obtaining a (hopefully) sparse vector $x$ we replace the matrix $A$ with $A-(x_1^TAx_1)x_1x_1^T$ and repeat the above steps to obtain the next sparse component values.

The question then is "when to stop?". Two approaches are proposed. First, at each iteration $i$, for all $i\lt j$, we include the constraint $x_i^TXx_i=0$ to make sure that each principal component we compute is orthogonal to the previous ones. Then the procedure stops after $n$ steps automatically (there will be no solution to the $n+1$'th problem).

The other way is stop as soon as all members of $A$ get less than $\rho$, because at that point elements of $A$ will be less than the noise value $\rho$.

### Solving the Semidefinite Problem

The cardinality constraint in the formulation (or its corresponding term in the penalized form) introduces a quadratic number of terms in the problem. This makes it practically impossible to use interior-point method to solve the problem for large values of the input dimension. So we need to use other existing methods for solving the problem, but then, there will be a matter of speed. Denoting the required accuracy by $\epsilon$, we can expect an interior-point-based program to converge after $\textstyle O(F(n)\log\frac1{\epsilon})$ iterations, for some function $F$ of input size $n$. Here we will manage to solve the problem using $\textstyle O(\frac{F(n)}{\epsilon})$ iterations using a first-order scheme, for some function $F$ of input size.

First-order method can solve a problem after $\textstyle O(\frac1{\epsilon})$ iterations, for some function $F$ of input size, if the problem satisfy a smoothness constraint. But in our case, $X\succeq 0$ is not smooth and as a result, an application of first-order method to the our problem will result in an algorithm stopping after $\textstyle O(\frac1{\epsilon^2})$ iterations, for some function $F$ of input size, which is too slow. To address this problem, we consider the formulation (5) and then we define a smooth approximation of the function $\lambda^{\max}$.

To come up with a smooth approximation of our goal function, we define $f_{\mu}(X)=\mu\log\textbf{Tr}(e^{\frac X\mu})$. Then, one can verify that $\lambda^{\max}(X)\leq f_\mu(X)\leq\lambda^{\max}(X)+\mu\log n$ and so, for $\textstyle\mu=\frac{\epsilon}{\log n}$, $f_{\mu}$ is a smooth approximation of $\lambda^{\max}$ with an additive error of $\epsilon$. This way we obtain a scheme for solving the program in $\textstyle\frac d{\epsilon}\sqrt{\log d}$ iterations, each taking $O(d^3)$ time.

## Experimental Results

Each point in figures below corresponds to an experiment on 500 genes. The points are pre-clustered to 4 clusters based some prior knowledge. The top three principal components are computed using each of PCA and Sparse PCA methods, and the points are plotted in the bases defined by these three components. For the PCA, each principal component is a combination of all 500 variables (corresponding to 500 genes) while in sparse PCA each involves variables corresponding to at most 6 genes.

Figure 1. Distribution of gene expression data in the PCA vs. Sparse PCA. The point colors are based on an pre-computed independent clustering.

In the next figure, the left diagram compares the cumulative number of non-zero elements in principal components in three methods: SPCA, the method we explaind with $k=5$, and the method we explaind with $k=6$. In the right diagram, the cumulative percentage of total variance explained by the first principle components resulted from PCA, SPCA (dashed) and the method we explained with $k=5$ and $k=6$ (solid lines).

Figure 2. Cumulative cardinality and total percentage of total variance explained by the first principle components resulted from PCA, SPCA (dashed) and the method we explained with $k=5$ and $k=6$ (solid lines).


text:
OpenGL - Bind Buffer Arrays question

So this is probably a simple question:

But how come sometimes I see it like this:

GLuint vbuffer;
glGenBuffers(1, &vbuffer);
glBindBuffer(GL_ARRAY_BUFFER, vbuffer);

but when done in multiples:

GLuint VBOs[2]
glGenBuffers(2, VBOs);
glBindBuffer(GL_ARRAY_BUFFER, VBOs[0]);

Im curious why glGenBuffers has a & in the first one, and not in the second case. Is it because arrays are always passed by reference? and is vbuffer a POINTER to something or just a simple data structure (GLint in this case).

• In C or C++, in this context, VBOs is short for &VBOs[0]. You could write &VBOs[0] if it confuses you. – user253751 Nov 5 '15 at 23:10

format: plain


text:
# Article

Full entry | PDF   (0.9 MB)
References:
[1] Kаlаs J.: Asymptotic behaviour of the solutions of the equation dz/dt =f(t, z) with acomplex-valued function f. Proceedings of tһe Internаtionаl Colloquium on Quаlitаtive Theory of Differentiаl Equаtions, August 1979, Szeged - Hungаry, Seriа Colloquiа Mаthemаticа Societаtis János Bolyаi & North-Hollаnd Publishing Compаny, pp. 431-462.
[2] Kаlаs J.: On the asymptotic behaviour of the equation dz/dt =f(t,z) with a complex-valued function f. Arch. Mаth. (Brno), 17 (1981), 11-22. MR 0672484
[3] Kаlаs J.: On certain asymptotic properties of the solutions of the equation ż = f(t, z) with a complex-valued function f. Czech. Mаth. Journаl, to аppeаr.
[4] Kаlаs J.: Asymptotic properties of the solutions of the equation ż = f(t, z) with a complex-valued function f. Arch. Mаth. (Bmo), 17 (1981) 11З-124.
[5] Ráb M.: The Riccati differential equation with complex-valued coefficients. Czech. Mаth. Journаl 20 (1970), 491-503. MR 0268452
[6] Ráb M.: Geometrical approach to the study of the Riccati differential equation with complex-valued coefficients. J. Differentiаl Equаtions 25 (1977), 108-114. MR 0492454
[7] Ráb M.: Asymptotic behaviour of the equation x" + p(t)x' + q(t)x = 0 with complex-valued coefficients. Arch. Mаth. (Bгno) 11 (1975), 193-204. MR 0404776

Partner of


text:
# PoS(LATTICE2016)060

Study of the sign problem in canonical approach

R. Fukuda, A. Nakamura, S. Oka, A. Suzuki, Y. Taniguchi

Contribution: pdf

Abstract

It is known that the sign problem emerges as a complex phase of the canonical partition function in canonical approach. We confirmed that the origin of the complex phase is breaking of the charge conjugation with a simple model. The main purpose of this study is suppressing the complex phase even at low temperature. We realized it by increasing the number of gauge field configurations. In particular, the complex phase was suppressed less than $\pi/2$ in $1\sigma$ for $N_B<7$ at $T=0.81T_c$


text:
# How to get PI based results in Octave?

In Octave

octave.exe:4> asin(1/2)
ans =  0.52360

is it possible to get the result in form of:

ans = pi/6

It seems not possible with format function, but may be there are any tricks.

-

I assume it's the same as with Matlab, where you would do sym(asin(1/2)).

-
This relies on the symbolic toolbox in matlab, so probably you will need a package like described here: math.stackexchange.com/questions/128812/… – Dennis Jaheruddin Nov 21 '12 at 12:31

An alternative might be to convert into units of $\pi$ explicitly, and then seek a rational approximation:

> rats(asin(1/2) / pi)
ans =       1/6
-

Try using "format rat". "rat" stands for rational, which puts answers in fractional terms. Using "format rat" puts ALL terms in fractional terms (unless you switch back to a different format such as "format g"), as opposed to rats([expression]), which (I assume - I've never used that) only puts [expression] in rational terms.

-


text:
# ld Eelectrohome DVD player MODEL EH8181A That stopped powering on!

Status
Not open for further replies.

#### Mr-Robby

##### New Member
Hi Folks,

I have a very old Eelectrohome DVD player MODEL EH8181A That stopped powering on today.

Thought I would post the Power Supply board to see if anyone knows what the most likely issue would be here.

I thought of things like:

1) Fuse - not sure how to identify where the fuse is here and how it could be replaced once located?
2) I see something that looks like burn on top of one of the what looks like a capacitor, could just be rust, if that blew, any way to replace it?
3) anything else to look at?

I can't seem to find my volt meter but I wouldn't know what to check even if I could.

Thanks! Rob

#### dknguyen

##### Well-Known Member
Is the cap bulging? It's hard to tell from the angle of the light in the photo but it looks different from the other caps which are for sure not bulging. If so, replace it.

#### Mr-Robby

##### New Member
Is the cap bulging? It's hard to tell from the angle of the light in the photo but it looks different from the other caps which are for sure not bulging. If so, replace it.
You mean the thing that looks like a capacitor and has discoloration on top? If so, yeah, I think it is bulging. Does that mean it blew up? How do I get a new one and how would I replace it? thanks!

#### dknguyen

##### Well-Known Member
You mean the thing that looks like a capacitor and has discoloration on top? If so, yeah, I think it is bulging. Does that mean it blew up? How do I get a new one and how would I replace it? thanks!
Just compare it with all the other caps next to it. There is normally a slight indentation or it's flat the most. If it's anything else it;s bulging. It hasn't blown up yet...but it's definitely not doing what it should be doing.

Get a soldering iron and heat up the cap's pins on the underside and it pull it out. Heat one pin at a time and tilt it out, or heat both at the same time if your tip is wide enough. You can refer to some youtube videos.

The capacitance and voltage rating should be written on the cap. Find another one of the same diameter and pin spacing with at least the same voltage and approximately the same capacitance. You can go to a local electronics shop if you have one and ask or order online. Mouser or Digikey are the main places I order form.

It's called an leaded electrolytic capacitor. If ordering online make sure to check datasheets for physical dimensions.

Last edited:

#### Mr-Robby

##### New Member
Ok, thanks. Will give it a try. I am new at playing aorund with these things, usually I just throw them in the trash and buy a new one but wanted to see if I could actual fix one for fun.

#### dknguyen

##### Well-Known Member
Ok, thanks. Will give it a try. I am new at playing aorund with these things, usually I just throw them in the trash and buy a new one but wanted to see if I could actual fix one for fun.
If it's really for fun, you could try powering it up after the cap is removed but before replacing it with a new one. It may run (depends how critical that capacitance was. It may have been working in tandem with the one right next to it and if the value wasn't critical it will still run, albight with degraded performance.

#### Ylli

##### Active Member
Electrolytic caps are polarized. Note any polarity markings on the old cap and be sure to install the new one with the same polarity.

#### Mr-Robby

##### New Member
Electrolytic caps are polarized. Note any polarity markings on the old cap and be sure to install the new one with the same polarity.
Ok, thanks V much for your insight.

#### Western

##### Member
Two more points ...

1) Most switched mode power supplies use capacitors with 105 degree rating ... not just the standard of 85 degrees ... so add that to the requirement for a replacement.

2) Again, with switched mode power supplies ... often when they will not operate ... the mains capacitor (large cylindrical cap - top rh corner of your photo) will often remain charged.

If you find your multimeter, measure for any voltage across it before you go too far ... it will help us in diagnosing more if replacing the damaged one doesn't solve your problem.

If it is still charged, it will pay to discharge it before you do to much else ... either using a resistor of say 5 watts and a resistance of a few hundred ohms ... or if you like sparks and loud noises ... just short it out with a screwdriver.

If you don't mind surprises and want reminding of the curse words you know ... just leave it charged.

#### Mr-Robby

##### New Member
Ok so I measured the little bugger best I could:

Diameter: 8mm
uF: 1000 uF
Temperature: 105 celcius
Height: 20mm
Pin Distance: 5mm
Voltage: 16V
Brand: <BH>
Polarity: TBD

So how would I ensure I get the right part and size?

Mouser or Digikey looks cool but I get lost on the first menu for filtering.

#### Mr-Robby

##### New Member
I took a stab at it, does this look right?

#### dknguyen

##### Well-Known Member
Ok so I measured the little bugger best I could:

Diameter: 8mm
uF: 1000 uF
Temperature: 105 celcius
Height: 20mm
Pin Distance: 5mm
Voltage: 16V
Brand: <BH>
Polarity: TBD

So how would I ensure I get the right part and size?

Mouser or Digikey looks cool but I get lost on the first menu for filtering.
Ok so I measured the little bugger best I could:

Diameter: 8mm
uF: 1000 uF
Temperature: 105 celcius
Height: 20mm
Pin Distance: 5mm
Voltage: 16V
Brand: <BH>
Polarity: TBD

So how would I ensure I get the right part and size?

Mouser or Digikey looks cool but I get lost on the first menu for filtering.
Does this help?
https://www.digikey.ca/products/en/capacitors/aluminum-electrolytic-capacitors/58?k=capacitor&k=&pkeyword=capacitor&pv508=4&FV=mu1000µF|2049,mu100V|2079,mu125V|2079,mu150V|2079,mu160V|2079,mu165V|2079,mu16V|2079,mu175V|2079,mu180V|2079,mu200V|2079,mu20V|2079,mu210V|2079,mu220V|2079,mu225V|2079,mu230V|2079,mu250V|2079,mu25V|2079,mu280V|2079,mu28V|2079,mu300V|2079,mu30V|2079,mu315V|2079,mu330V|2079,mu350V|2079,mu35V|2079,mu360V|2079,mu385V|2079,mu400V|2079,mu40V|2079,mu415V|2079,mu420V|2079,mu42V|2079,mu450V|2079,mu45V|2079,mu475V|2079,mu500V|2079,mu50V|2079,mu525V|2079,mu550V|2079,mu55V|2079,mu56V|2079,mu575V|2079,mu580V|2079,mu600V|2079,mu60V|2079,mu630V|2079,mu63V|2079,mu650V|2079,mu70V|2079,mu71V|2079,mu75V|2079,mu80V|2079,ffe0003a,1c0002,1c0006,402e23,1140050,1f140000&quantity=0&ColumnSort=0&page=1&pageSize=100

BTW, Mouse prices are cheaper but the Digikey search is superior.

You need not enter all search criteria in at the same time. Enter the most important things and filter. It will remove narrow down and remove search criteria that aren't available for products based on the previous search criteria you already used.

#### dknguyen

##### Well-Known Member
I took a stab at it, does this look right?

View attachment 113087
The third one seems to have the proper lead spacing. The first one does not and the second its not listed for the second one in the search database (which means you have to check the datasheet but capacitor datasheets are overwhelming if you're starting out since a single datasheet covers multiple sizes and models, all with slight differences and only discernible by a vague code. The simpler the component, the more complicated the datasheet is the way it usually goes.).

If you use the search enough you learn not to always trust the search since sometimes information is entered in or are classified incorrectly, or not classified so they dissapear when filters are applied when they should appear...but let's not worry about that for now.

#### Mr-Robby

##### New Member
The third one seems to have the proper lead spacing. The first one does not and the second its not listed for the second one in the search database (which means you have to check the datasheet but capacitor datasheets are overwhelming if you're starting out since they cover every model available in that line)
yeah, i forgot to look at spacing column. So looks like #3 wins.

#### dknguyen

##### Well-Known Member
Looks good to me. Might want to throw some extra stuff like a few tools into the order because $1 capacitor for$8 worth of shipping through the mail seems like a waste.

Maybe solder, an iron, and hemostats if you don't already have them.

#### Mr-Robby

##### New Member
Looks good to me. Might want to throw some extra stuff like a few tools into the order because $1 capacitor for$8 worth of shipping through the mail seems like a waste.

Maybe solder, an iron, and hemostats if you don't already have them.
I plan to take these specs across to the Chinese mall here. There's a guy there with a bunch of old electronics and he may be able to easily source this part. will see. thanks everyone!


text:
# Probability that a vector in $\mathbb{Z}^n$ is primitive

A vector $v \in \mathbb{Z}^n$ is primitive if there does not exist some vector $v' \in \mathbb{Z}^n$ and some $k \in \mathbb{Z}$ such that $v = k v'$ and $k \geq 2$.

For a paper I'm writing right now, I'd like to know that a "random" vector in $\mathbb{Z}^n$ is primitive. Let me make this precise.

Let $\|\cdot\|_{1}$ be the $L^{1}$ norm on $\mathbb{Z}^n$, so $\|v\|_1 = \sum_{i=1}^n |v_i|$, where the $v_i$ are the components of $v$. Define $\mathcal{V}_k$ to be the number of vectors $v$ in $\mathbb{Z}^n$ such that $\|v\|_1 \leq k$. Define $\mathcal{P}_k$ to be the number of primitive vectors $v$ in $\mathbb{Z}^n$ such that $\|v\|_1 \leq k$.

I then want $$\lim_{k \rightarrow \infty} \frac{\mathcal{P}_k}{\mathcal{V}_k} = 1.$$ Assuming this is true, is there any nice estimate as to how fast it approaches $1$?

-
Even if you didn't know the exact answer to this question, the set of all vectors whose elements are even is a subset of non-primitive vectors of density 1/2^k, so the answer is manifestly less than 1. By inclusion-exclusion, it's at most (1 - 1/2^k)(1 - 1/3^k)(1 - 1/5^k)... for any finite collection of primes and at this point you might conjecture the relation with the zeta function. –  Qiaochu Yuan Nov 20 '10 at 12:29
@Qiaochu: Yes, that's the intuitive way to think about this. The probability of the vector not being a multiple of p is $1-1/p^n$. Inclusion exclusion or independence (from CRT) gives $\prod_p(1-1/p^n)$ for the probability of not being a multiple of any prime from a finite set. This quickly gives $1/\zeta(n)$ as an upper bound, but doesn't immediately give the lower bound or a rate of convergence. Expanding the product gives the expression in terms of the Mobius function in my answer from which, with a bit of work, you get the rate of convergence. –  George Lowther Nov 20 '10 at 13:04
@George: right. Using the argument given, for example, at qchu.wordpress.com/2010/11/09/… , you can show that if the natural density exists it must be 1/zeta(n), but you certainly need to do extra work to actually get that the density exists. –  Qiaochu Yuan Nov 20 '10 at 13:30

In the case where $n = 2$, you're asking for the "probability" that two integers are relatively prime; this is well-known to be $6/\pi^2$, not 1. In the general-$n$ case, the probability that $n$ integers are relatively prime is $1/\zeta(n)$.

Reference: http://en.wikipedia.org/wiki/Coprime

-
Thanks! Do you have a reference for this? –  T_P Nov 19 '10 at 20:20
@T_P : I have added a reference –  anonymous Nov 19 '10 at 20:23
Thanks, Chandru. I couldn't add a reference because I was teaching! –  Michael Lugo Nov 19 '10 at 21:07
No, problem sir! –  anonymous Nov 19 '10 at 21:29

Further to Michael's answer, not only does $\mathcal{P}_k/\mathcal{V}_k\to1/\zeta(n)$, but we can calculate a bound for the rate of convergence. I'll also give an argument which is a little different from the one given in his Wikipedia reference.

Noting that the set $\lbrace v\in\mathbb{R}^n\colon\Vert v\Vert_1\le k\rbrace$ has volume $ck^n$ (for a constant c depending only on the dimension n) and surface area proportional to $k^{n-1}$ gives $$\mathcal{V}_{k}-1=ck^n + O(k^{n-1}).\qquad\qquad{\rm(1)}$$ The '-1' on the left hand side is not relevant for large k as it can be absorbed into the O(kn−1) error term, and is just there so that (1) is also valid for small k < 1. Noting that every nonzero $v\in\mathbb{Z}^n$ decomposes uniquely as $v=mv^\prime$ for integer $m\ge1$ and primitive $v^\prime\in\mathbb{Z}^n$ leads to the following relation between $\mathcal{P}_k$ and $\mathcal{V}_k$, $$\mathcal{V}_k-1=\sum_{m=1}^\infty\mathcal{P}_{\frac{k}{m}}.$$ This can be inverted via the Möbius function μ, $$\mathcal{P}_k=\sum_{m=1}^\infty\mu(m)(\mathcal{V}_{\frac{k}{m}}-1).$$ In dimension $n > 2$, substituting (1) into this expression gives $$\mathcal{P}_k=\sum_{m=1}^\infty \mu(m)c k^n m^{-n} + O(k^{n-1}).\qquad\qquad{(2)}$$ The $O(k^{n-1})$ comes from the sum $\sum_m (k/m)^{n-1}$ from the remainder term of (1) which, for $n > 2$, gives $k^{n-1}$ multiplied by a convergent sum. Dividing through by $\mathcal{V}_k$, $$\mathcal{P}_k/\mathcal{V}_k=\sum_{m=1}^\infty\mu(m)m^{-n}+O(1/k)=1/\zeta(n)+O(1/k).$$ Edit: The case for $n=2$ is actually a little bit different, and we do not obtain such a good convergence rate. As the sum $\sum_m(k/m)^{n-1}$ does not converge, the error term in (2) does not apply. Instead, we can use $O(1_{\lbrace k\ge1\rbrace}k+1_{\lbrace k < 1\rbrace}k^2)$ for the error term in (1). This leads to an error of order $k\sum_{m\le k}m^{-1}+k^2\sum_{m > k}m^{-2}\sim k\log k$ in (2), giving $$\mathcal{P}_k/\mathcal{V}_k=1/\zeta(2)+O(\log k/k).$$ You can also look at the paper On the probability that k positive integers are relatively prime.

-


text:
# 2 ways to solve an ODE, different solutions…

Given the ODE $$y'= \left( \frac{x+2y+1}{2x-6} \right)$$

I have tried two methods to solve this equation. One method is seperate variables: $$(2x-6 )dy = (x+2y+1)dx$$ integrating gives $$2(xy-3y)=x^2+2xy+x+c$$, simplifying we get $$y= -(x^2+x+c)/6$$\

But on the other hand, if I try to solve by using an integrating factor u(x), I get a different solution. The method I used was rewrite in standard linear form, find integrating factor by calculating $u(x) = e^(\int p(x)dx)=1/(3-x)$. In the end I found that the solution was $y=(3-x)(2/(x-3)+1/ log(x-3) + C$

Whats going on?

-
Did you do a sanity check and take the derivative of both answers to see if either is correct? –  tpg2114 Nov 23 '12 at 12:20
@Hempo: You can't use the first method! what is $\int ydx$? –  Dennis Gulko Nov 23 '12 at 12:23
Separation of variables applicable here: $Y(y)\mathrm d\,y=X(x)\mathrm d\,x$. –  FrenzY DT. Nov 23 '12 at 13:05
Ok. The correct reasoning is that y(x) is a function of x, and hence it cannot be seen as a constant if you integrate it –  MSKfdaswplwq Nov 25 '12 at 13:10
Remember that $y$ is a function of $x$, so $\int ydx\neq xy$ and $\int xdy\neq xy$.


text:
# Rule with if and average value

in the rule I have this if:

``````val avg = ModbusPVPower850_ValueasNumber.averageSince(now.minusMinutes(5))

if (avg < 2500) {
Relay_2_status.postUpdate(avg.toString)
}
``````

Problem is, that this isn’t working
if I will do only this without if:

``````
Relay_2_status.postUpdate(avg.toString)

``````

it’s working, so problem must be in if declaration, but don’t know which

I suppose you could use `logInfo()` to find out what `avg` is - might there be a unit as part of it?

ok, let me chek later, now I don’t have a access to the terminal.

ok, I found it :

``````val PV_average = ModbusPVPower850_ValueasNumber.averageSince(now.minusMinutes(5)) as Number
``````

now it works


text:
## 2011-06-26

### Physicists are not physicians

I've just realized to have confused often the word "physicist" with the word "physician". Luckly, after a fast check, I was not so bold to label myself as a "physicist" (since I've currently discontinued physics studies alas) in cv or sites' profiles, otherwise I would have written "physician" instead. Though, in forums and blogs here or there maybe it's happened. So forgive me, it is very likely that where you read the word "physician", I meant "physicist".

## 2011-06-23

### Who gets it right? Enterprise encoding hell

One of my preferred arguments about computers is character encodings. Years ago I began to write several lectures about computers and programming, and the first were all about this subtle issue.

The fact is that people fail to understand that computers do not understand, they just execute codes (that is, data which can be interpreted as code by a special interpreter, the processor); they are just programmed to interpret data in a particular way and maybe to represent them some way. But there's no meaning in data themselves. So if you get the wrong interpreter for some data, what you get is garbage, or a polite error message. If you have an interpreter but you feed it with data it can't interpret, you get again garbage as output, or a polite error message.

What we read as letters for a computer are data (bytes), no more no less. Data that can be interpreted as letters, but in order to do so, the interpreter must know (be programmed to use) a sort of map between the data and the symbols shown on the screen, on the mobile phone display, or wherever.

Nowadays the ASCII "standard" is so widespread that we can disregard any other standard (like EBCDIC) though likely it's still alive somewhere. The ASCII encoding "states" that the byte XY, when interpreted as a "letter", must be shown as the letter "L"; put numbers (another interpretation of bytes) in stead of XY and other letters in stead of L and you have a way of interpreting a sequence of bytes as a sequence of letters. Indeed ASCII has also bytes that have no "printable" symbols associated; e.g. there's a byte that means go to the next line (new line), or put a space (space) and so on.

Anyway, ASCII has no too much special symbols, accented letters (since it was born as an american/english stuff only) and so on. So the pain starts when you want to see that nice symbol on every computer interpreting your data as a document made of readable symbols.

The encoding babel begins at that point. It could be interesting to trace the history of character encodings. Anyway for any practical purpose it is more interesting to aknowledge the fact that there were great efforts to reduce this babel and in the same time to increase the number of symbols representable in a "unified" manner.

But there's no way to escape from the choice of the correct data to be fed to an interpreter that expects those data.

So, if you as programmer, write in your preferred IDE or editor, a line like this:

printf("il costo è di 10€"); //avoid this!

or if you feed your DB with a statement that looks to you like

INSERT INTO aTable VALUES('il costo è di 10€', 10); -- avoid this!

and you know nothing about the encoding your IDE/editor uses and your DB accepts/translates into, and the encodings that are understood by the systems where you suppose to display those texts, then you are getting it badly wrong, to say the least.

About DB, when you want to be sure it does not "scramble" your text, you must force it to believe it is data, not text. You loose something (ordering capabilities according to specific collating sequences), but you are sure that the bytes you wrote are "transmitted" as they were to the final end-point. (The question about in which encoding those bytes represent "è" and "€" remains, i.e. you must input those symbols in a desidered known encoding, or input them as bytes, if possible).

Now, some people may ask how it is possible that I wront "è€" and the reader reads "è€"? HTML (which is ASCII based) has special "metainformartions" that are used to specify the encoding of the content. If the encodings match, I succeeded showing you the symbols I see: we are saying to the "interpreter" (the browser that shall choose the right symbol according to the input) which "map" to use.

The default encoding for web pages if nothing is specified is ISO-8859-1, AKA Latin1, which is 8bit superset of the ASCII encoding and a single byte encoding — the first 128 codes of Latin1 are the same of the ASCII (which has "only" 128 codes). (Instead of ISO-8859-1, it could be ISO-8859-15 aka Latin15; this has only the "currency symbol" replaced by the "euro symbol").

Anyway, the most general encoding is UTF(8) which is a way to encode the biggest "map" available, the "map" that tries to hold every symbols humans have created (UCS i.e. Universal Character Set). Modern O.S. uses this encoding as default encoding, but of course it is not always true.

To get a little bit deeper in the matter, wikipedia has articles about the subject. It is worth reading a lot about the topic even by employees of big enterprises.


text:
By using the substitution $p=x+\frac{1}{x}$, show that the equation $$2x^4+x^3-6x^2+x+2=0$$ reduces to $2p^2+p-10=0$.

I can't think of anything that produces a useful result, I tried writing p as $p=\frac{x^2+1}{x}$ and finding areas to substitute but have come with no progress. Could someone offer a slight hint on how to proceed?

Thanks

-

Like this,

as the equation is Reciprocal Equation of the First type with $x\ne0,$

divide either sides by $\displaystyle x^{\frac42}=x^2$ to reduce the degree of the eqaution by half

$$2x^2+x-6+\frac1x+\frac2{x^2}=0$$

$$\implies 2\left(x^2+\frac1{x^2}\right)+\left(x+\frac1x\right)-6=0$$

$$\implies 2\{\left(x+\frac1x\right)^2-2\}+\left(x+\frac1x\right)-6=0$$

Reference : Reciprocal Equation is explained here:

Chapter XI of Higher Algebra,Barnard & Child and

Article $568−570$ of Higher algebra, Hall & Knight

-
how do I know how to write $x^2+\frac{1}{x^2}$ as $(x+\frac{1}{x})^2-2$? –  salman Aug 30 '13 at 18:03
@user90771, as you have $p=x+\frac1x,$ and we know $(x+\frac1x)^2=x^2+\frac1{x^2}+2$ –  lab bhattacharjee Aug 30 '13 at 18:07
also jee, how do I know that I must divide by $x^2$? thanks –  salman Aug 30 '13 at 18:08
@user90771, please find the edited answer, hope this should explain the idea –  lab bhattacharjee Aug 30 '13 at 18:11
this is an amazing resource thanks bhai –  salman Aug 30 '13 at 18:24

Hint: Try working backwards. Start with: $$2\left(x + \dfrac{1}{x}\right)^2 + \left(x + \dfrac{1}{x}\right) - 10 = 0$$ then try to obtain the original equation by expanding then clearing the denominators.

-

$$p=x+x^{-1}$$ $$p^2=x^2+2+x^{-2}$$ $$2x^4+x^3-6x^2+x+2=0$$ $$2x^2+x-6+x^{-1}+2x^{-2}=0$$ $$2(x^2+2+x^{-2})+(x+x^{-1})-10=0$$ $$2p^2+p-10=0$$

-
what indicates that I must divide by $x^2$? –  salman Aug 30 '13 at 18:06
@user90771 It's the standard usual method which works well. –  namehere Aug 30 '13 at 18:07
@user90771 You divide by x square because it just happens to work. In general quartic equations are not subject to substitutions as easy as this one. –  namehere Aug 30 '13 at 18:08

\begin{align} 2x^4+x^3-6x^2+x+2 & = x^2\Big(2x^2 + x - 6 + \frac1x + \frac{2}{x^2}\Big) \\[12pt] & = x^2 \Big(2\left(x+\frac1x\right)^2 + \left(x+\frac1x\right) - 10 \Big) \\[12pt] & = x^2(2p^2 +p-10). \end{align} This equals $0$ only if $x=0$ or the second factor equals $0$. But clearly $x=0$ is not one of the solutions.

-


text:
# Calculating a Positive Relative error

I need to graph the error between a set of observed values and an expected values.

I calculate the relative error for each point as:

(observed - expect) / expected

This error however bounces up and down because the behaviour is sinusoidal. I want to make this error always a positive quantity. Is there an error measure like this? if so, what's its name?

I have looked for error measures, but all I find is absolute and relative, and I am interested in a positive relative error.

Any help appreciated, Ted.

-
you can always use the absolute value. –  picakhu Feb 2 '11 at 20:27
for example, if observed-expected < 0, then use expected-observed instead –  picakhu Feb 2 '11 at 20:28
Does this have a name? –  Ted Flethuseo Feb 2 '11 at 20:38
@Ted Flethuseo: Yes, absolute value, as mentioned by picakhu. It's written $|x|$, defined by $$|x|=\left\{\begin{array}{ll}x &\mbox{if x\geq 0}\\-x &\mbox{if x\lt 0.}\end{array}\right.$$ Most computers/spreadsheets can compute it, usually with the abs command. –  Arturo Magidin Feb 2 '11 at 20:41
absolute value.. –  picakhu Feb 2 '11 at 20:41

$$\frac{(O - E)^2}{E}$$
This is related to the test statistic for Pearson's $\chi^2$ test.


text:
## The Curse of Cause of Death Models

Stephen's earlier blog explained the origin of the very useful result relating the life-table survival probability $${}_tp_x$$ and the hazard rate $$\mu_{x+t}$$, namely:

${}_tp_x = \exp \left( - \int_0^t \mu_{x+s} \, ds \right). \qquad (1)$

To complete the picture, we add the assumption that the future lifetime of a person now aged $$x$$ is a random variable, denoted by $$T_x$$, and the connection with expression (1) which is:

${}_tp_x = \Pr[ T_x > t ]. \qquad (2)$

The 'package' of random lifetime $$T_x$$, hazard rate $$\mu_{x+t}$$, survival function $${}_tp_x$$, and expression (1) tying everything together, sums up the mathematics of a survival model.  For the statistics, we have observations on the $$i^{\rm th}$$ of $$n$$ individuals, namely $$t_x^i$$ the observed lifetime, and $$d_x^i$$ an indicator equal to 1 if observation ended in death and 0 otherwise.  We suppose each $$(t_x^i,d_x^i)$$ is a sample drawn from the distribution of a bivariate random variable $$(T_x,D_x)$$.  The target of estimation is the distribution of $$(T_x,D_x)$$.

What does the 'package' above look like if we try to model cause of death, as well as time of death?  A natural extension is to associate a random future lifetime with each separate cause.  As a simple example, suppose heart attack and stroke are the only causes of death.  Then, for a person aged $$x$$, define:

\begin{eqnarray*} U_x & = & {\rm the\ time\ until\ death\ by\ heart\ attack\ occurs} \\ V_x & = & \mbox{the time until death by stroke occurs}. \end{eqnarray*}

The two causes of death compete with each other, to be the first to claim their victim, hence the name competing-risks model for models of this kind.  We discuss these in more depth in Chapter 16 of our forthcoming book, Modelling Mortality with Actuarial Applications; Crowder (2001) is another good source.

We can never observe both $$U_x$$ and $$V_x$$, but only their minimum, which we denote by $$T_x^*$$:

$T_x^* = \min[U_x , V_x].$

Associated with $$T_x^*$$ is a survival function denoted by $${}_tp_x^*$$:

${}_tp_x^* = \Pr[ T_x^* > t ].$

Define a hazard rate $$\lambda_{x+t}$$ associated with $$U_x$$ (in the presence of $$V_x$$) and a hazard rate $$\nu_{x+t}$$ associated with $$V_x$$ (in the presence of $$U_x$$) in exactly the same way as the hazard rate $$\mu_{x+t}$$ is associated with $$T_x$$ (see Stephen's earlier blog). Then we can show that:

${}_tp_x^* = \exp \left( - \int_0^t ( \lambda_{x+s} + \nu_{x+s} ) \, ds \right). \qquad \mbox{(3)}$

So we have the mathematical 'package' of random lifetimes $$U_x$$ and $$V_x$$ and their observable minimum $$T_x^*$$, hazard rates $$\lambda_{x+t}$$ and $$\nu_{x+t}$$, survival function $${}_tp_x^*$$, and expression (3) tying everything together.

The statistical side of the competing risks model is much the same as before, because we can observe only $$T_x^* = \min[U_x,V_x]$$, and the cause of death (or, at least, what prevailing practices define as the cause of death). Our observations of the $$i^{\rm th}$$ of $$n$$ individuals will take the form of the time of death $$t_x^{i*}$$, and an indicator $$d_x^{i*}$$ of cause of death, for example $$d_x^{i*} = 1$$ if heart attack caused death or $$d_x^{i*} = 2$$ if stroke caused death.  The target of estimation is the joint distribution of the bivariate random variable $$(T_x^*, D_x^*)$$ from which the observed data $$(t_x^{i*}, d_x^{i*})$$ are supposed to be sampled.

Now comes the crux — we have said nothing about the dependence or independence of $$U_x$$ and $$V_x$$.  For given hazard functions $$\lambda_{x+t}$$ and $$\nu_{x+t}$$, the lifetimes $$U_x$$ and $$V_x$$ may be independent, or may have any kind of dependence that is mathematically possible.  Expression (3) is still true, regardless.  But this means that the distribution of $$(T_x^*, D_x^*)$$ is always the same, whatever the dependence between $$U_x$$ and $$V_x$$ may be.  The time to death is governed by the same survival function $${}_tp_x^*$$, and if death occurs at time $$t$$ the probability that the cause was a heart attack is always the same, namely:

$\Pr[\, \mbox{Heart attack} \mid \mbox{Death at time } t \,] = \frac{\lambda_{x+t}}{\lambda_{x+t} + \nu_{x+t}}.$

As a result, it is impossible to tell, from any amount of data, whether these two competing causes of death are dependent or independent. And, if those causes of death are dependent, it is equally impossible to say what the nature of that dependence is.

Our chosen causes, heart attack and stroke, very likely are strongly related, because of genetic and lifestyle factors (e.g. smoking and alcohol consumption).  We would expect a cause-of-death model to let us estimate such dependencies.  But no.

This is the identifiability problem (or, as one statistician called it, crisis).  It is usually attributed to Tsiatis (1975), although it was noticed earlier.  It poses an obvious challenge to cause-of-death modelling of a rather fundamental nature.  Dependencies between different causes, which we believe to exist on biological grounds, can never be estimated in a survival model based on random lifetimes.  They are forever out of reach, unless we equip the model with an additional mechanism giving rise to the dependencies, which often amounts to assuming the very thing we set out to estimate.  A curse indeed.

References

Crowder, M. (2001). Classical Competing Risks. Chapman & Hall CRC, Boca Raton, FL.

Macdonald, A. S., Richards. S. J. and Currie, I. D. Modelling Mortality with Actuarial Applications. Cambridge University Press (forthcoming).

Tsiatis, A. A. (1975). A Nonidentifiability Aspect of the Problem of Competing Risks. Proceedings of the National Academy of Sciences, U.S.A., 72, 20–22.